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ABSTRACT 


This thesis applies time-frequency transformations to radar signals. Specifically, it 
considers the feasibility of applying time-frequency transformations to extract the intra¬ 
pulse modulation parameters of radar signals. In this work, we consider radar signals with 
analog pulse compression; specifically linear or hyperbolic intra-pulse modulation. 
Several time-frequency transformations are investigated to identify which one gives the 
most accurate image representation for signals in noisy environments. Next, image 
processing techniques are applied in conjunction with an adaptive curve fitting method, 
for the hyperbolic modulation scheme, to extract the parameters of the frequency 
equation. Results show that for the linear chip case the frequency equation can be 
estimated with small error down to SNR equal to — lOdB. The proposed method for the 
hyperbolic chirp modulation is less immune to noise degradation and it can be used down 
to SNR level equal to 2dB. 
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I. INTRODUCTION 


The analysis of a signal in both the time and frequency domain gives a more 
complete description than that provided by individual analysis in either domain. This 
advantage was recognized a long time ago but the complexity of the joint time-frequency 
algorithms made their implementations an extremely slow operation due to computing 
power limitations. However, as computers become more powerful, joint time-frequency 
analysis may be applied in a larger number of real-world problems. 

A. OBJECTIVE 

This thesis investigates the application of time-frequency transformations to the 
extraction of intra-pulse modulation parameters from radar signals. This information can 
be a useful tool to identify the specific radar type and help with jamming techniques, if so 
desired. Intra-pulse modulation parameters are usually extracted using hardware schemes 
which are well suited to extract the instantaneous frequency at high SNR levels. 

Our study considers the problem from a different angle and investigates the 
application of time-frequency and a basic image processing technique to extract the 
information. Joint time-frequency transformations can be implemented inexpensively and 
give accurate results in medium to high SNR levels. We consider radar signals with 
analog pulse compression; specifically linear or hyperbolic intra-pulse modulation. 
Several time-frequency transformations are investigated to identify that which gives the 
most accurate image representation for signals in noisy environments. Next, image 
probessing techniques are applied used in conjunction with an adaptive curve fitting 
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method, to extract the parameters of the frequency equation. The type of modulation, as 
well as, the start and stop time of the pulse under investigation was assumed to be known. 

B. THESIS ORGANIZATION 

Radar signals with large time bandwidth products are introduced in chapter n. 
Chapter m presents an overview of the time frequency methods considered. Chapter IV 
discusses the basic idea behind the Radon transform and its application in detecting line 
parameters from a noisy image. Chapters V and VI present the complete methods and the 
simulations conducted to test the schemes derived to extract the intra-pulse modulation 
parameters for linear and hyperbolic modulation. Finally, conclusions and 
recommendations for further research are presented in chapter VH. 
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II. RADAR SIGNALS WITH LARGE TIME BANWIDTH PRODUCTS 


Ideally we would like a radar to produce good range resolution as well as large 
detection range. Unfortunately, these requirements are contradictory since the first 
requires a small pulse period, while the second requires a large pulse duration. The 
problem can be solved with a technique that uses intra-pulse modulation or, as it is better 
known, pulse compression. This chapter present the basic ideas behind pulse 
compression. 


A. RESOLUTION 

A very important aspect of the radar is its resolution properties. By this we mean 
the ability to detect multiple targets that are close together. Generally, there are four types 
of resolution: 

- Range resolution 

- Horizontal (azimuth) cross-range resolution 

- Vertical (elevation) cross-range resolution 

- Doppler frequency resolution. 




Figure 1: Range and Cross-Range Resolution From [1] 
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1. Range Resolution 

Range resolution is the ability to distinguish between targets, which are at the 
same angular position but are at different distances to the radar. This resolution is directly 
proportional to the radar signal bandwidth: the larger the bandwidth is, the better the 
range becomes. In order to resolve two targets, the basic criterion is that they must be 
separated by at least the range equivalent of the width of the processed echo pulse as 
illustrated in Figure 2. 



Processed Unresolved Near resolved Resolved 

pulse width 


a. Wide processed pulse width 



Processed Unresolved Near resolved Resolved 

pulse width 


Figure 2: Range Resolution and Processed Pulse Width. From Ref. [1 ] 

2. Horizontal and Vertical Cross-Range Resolution. 

Cross-range resolution is defined as the ability to distinguish between targets that 
have the same distance from the radar but are at a different azimuth (horizontal) and/or 
different elevation (vertical). This resolution is inversely proportional to the antenna 
beamwidth. That means the smaller the beamwidth is, the better the cross-range 
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resolution is. The basic constraint is that the targets must be separated by at least the 


antenna beamwidth. Figure 3 shows the effect of the antenna’s beamwidth to the 
horizontal cross-range resolution. 


Antenna 



Targets separated by 
less than beamwidth 



Targets separated by 
antenna beamwidth 


Targets separated by 
greater than beamwidth 


Signal echoes - 
Unresolved 


Signal echoes - Signal echoes - 

Near resolved Resolved 


Figure 3: Antenna Beamwidth and Horizontal Cross-Range Resolution. From 
Ref. [1] 


3. Doppler Resolution 

Radar doppler resolution is defined as the ability to distinguish between targets 
that have the same range, azimuth and elevation but have different radial velocities. The 
basic criterion is that the Doppler frequencies of the targets must be separated by at least 
one cycle over the time of observation. 


B. PULSE COMPRESSION 

A radar should have a large a detection range. If we assume that we use the best 
techniques for antenna and signal processing gain, the only way to increase the range is 
by increasing the energy of the transmitted signal. There are two ways to increase the 
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energy; by increasing its amplitude (power) of the pulse, or by increasing the pulse 
duration. Since all transmitters have peak power limitations, the most common choice is 
to increase the pulse duration, which is in contradiction to the need for narrow pulses to 
ensure good range resolution. The technique that uses long pulses with high bandwidth is 
called pulse compression. 

1. Analog Pulse Compression 

This technique uses modulated radio frequency signals. Let’s assume that we have 
a finite pulse duration radar signal with a constant radio frequency defined as: 

c(t) = r(t)-s(t), (H-l) 

where r(t) is a pulse and s(t) is a sinusoid function. Let C(f), R(f) and S(f) be the 
Fourier transform of the time domain signals, then: 

= (H-2) 

where the operator (*) denotes the convolution operation. Thus, the spectrum of 
the total signal has the form of a sine function centered at the frequency of the sinusoid 
with a null to null bandwidth equal to the double of the reciprocal of the width of the 
pulse in the time domain, as illustrated in Figure 4. 
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Figure 4: Bandwidth of a CW pulse. 


The Fourier transform S(f) is not a delta function when the radio signal s(t) is a 
modulated sinusoid. This in turns affects the bandwidth of the transmitted signal and 
gives a null to null bandwidth greater than 2/T, as shown in Figure 5. 



Figure 5: Bandwidth of a linear modulated pulse. 
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Next we will examine the two most widely used types of pulse intra-modulation; 
linear and hyperbolic modulations. 

a) Linear Modulation 

In this case the frequency of the radio signal inside the pulse changes 
linearly with a slope k from an initial frequency f 0 which leads to: 

/(0 = /o +k-t. (n-3) 

Where/o is the initial frequency at t=0.Thus, the transmitted signal c(t) is 

given as: 

c(t) = rect (~) • exp( j2aj fdt) =rect (p) • ^' 2 *(/o' + f' 2 > . (H-4) 

b) Hyperbolic Modulation 

In this case, the signal inside the pulse is modulated following a hyperbolic 
equation. The resulting frequency equation is given as: 

f(t) = j + b. (H-5) 

Note that the above time varying frequency f(t) is not finite for t=0. Thus, 
the above definition is modified by introducing a time shift t 0 to ensure a finite frequency 
value for t=0, which leads to: 

= + (n-6) 

* +r 0 
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Thus, the instantaneous frequency ./ft) can be viewed as being a part of the 


hyperbola defined from the equation (II-5) with an initial value of f 0 = 




Time 

Figure 6: Hyperbolic frequency modulation. 


As a result, the equation of the transmitted signal is given by: 

c(t) = rect (j) • exp( j2n\ fdt) =rect (j) ■ . (n-7) 


2. Digital Pulse Compression 

Digital pulse compression involves phase-coded waveforms. These are usually bi¬ 
phase modulated sinusoids with the two possible phases to be at 0° and 180°. The overall 
digital waveform consists of an array of N subpulses, each one with an initial phase of 0° 
or 180° and a time length x s . The width of the total waveform is Xe- Figure 7 shows a 
biphase coded signal using a 13-bit biphase sequence. In this figure a “+” corresponds to 
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a subpulse with an initial phase 0° and a corresponds to a subpulse with an initial 
phase 180°. 



a. Waveform 



± 

+ + + 

0 

□ 

□ 

[+1+1 

- 4- - 

m 


- 

T s 

U- 

t e 

-- — -- 



Figure 7: Phase coded waveform. From Reffl] 

The codes that are used in the phase coded waveforms should have a delta-like 
autocorrelation function, to allow pulse compression to take place without windowing the 
signal. The two most widely used codes are the Barker codes and the Pseudorandom 
codes [1]. 
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III. TIME-FREQUENCY REPRESENTATIONS 

Generally, a signal can be represented using different types of decompositions. 
The two most widely used representations are the time domain and the frequency domain 
representation. The first shows how the amplitude of the signal changes with respect to 
time while the second shows how often these changes occur. These representations are 
uniquely related to each other via the Fourier transform. 

A. FOURIER ANALYSIS 

1. Fourier Series 

The basic idea behind the Fourier decomposition is to express a periodic signal as 
a weighted sum of sinusoids. Let’s assume that we have a signal x(t) with period T. 
Provided that it satisfies the Dirichlet’s conditions, then it can represented as follows: 

x(t)='Y J a k -e M2nm \ (m-i) 

£=-~ 

where 

. a k =— j x(t)-e~ jk(27rlT)t dt .* (IH-2) 

The frequency 1/T is called the fundamental frequency. The equation (UI-1) gives 
the general case where the type of basis functions are complex exponentials. 

2. Fourier Transform 

By analogy with equation (ID-2), the Fourier transform X(f) of a continuous signal 
x(t) is given by: 


11 



(ffl-3) 


*(/) = jx(t)-e~ j 2 ) *dt. 

— oo 

The plot of the magnitude squared of the Fourier transform is proportional to the 
spectrum of the signal x(t). The original signal can be recovered from its Fourier 
transform with the inverse Fourier transform: 

oo 

x(t)= jX(f)-e~ J2 ^df . (HI-4) 

— oo 

The fast discrete version of Fourier transform revolutionized the signal processing 
area. However, there exist many signals which cannot be described accurately with this 
transformation. Although the spectrum indicates the frequencies and intensities of the 
signal, it provides no information regarding the variations of the signal characteristics as a 
function of time. Thus, the Fourier transform does not analyze non stationary signals very 
accurately. For example, consider two signals Xj(t) and x 2 (t), observed during the same 
time period T, as shown in Figure 8. The first one is the sum of two pure tones while the 
second one consists of the first tone for the first half time and then switches to the second 
tone for the rest of its duration. Their Fourier transform magnitudes look almost identical, 
however they are not exactly the same since the Fourier transform is a reversible 
operation. Figure 8 shows that it is not possible to extract specific information regarding 
the frequency variations as a function of time. Unfortunately, most signals in real life are 
not stationary, such as for example speech, music, vibration signals and many more. 
Analyzing non-stationary signals requires a tool which allows to represent the variations 
of the frequency content as a function of time; i.e., a joint time-frequency representation. 
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Several techniques have been developed for this type of representation. In this thesis, we 
will deal with two major classes: atomic decompositions and energy distributions. 


x1(t) 



Figure 8: Fourier transform of two signals. 

B. ATOMIC DECOMPOSITIONS 

The Fourier transform can be considered as the projection of a signal into an 
infinite set of sinusoids. The problem with sinusoids is that they are not localized in time, 
although they are perfectly localized in frequency. The basic idea behind the atomic 
decomposition is to decompose a signal into a set of atoms (i.e. functions) that are well 
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localized both in time and in frequency. These types of decompositions are called atomic 
decompositions or linear time-frequency representations. 

1. Short-Time Fourier Transform 

This method is an extension of the Fourier transform, where localization in time is 
introduced by windowing the sinusoidal functions used in the decomposition. Thus, the 
Sort-Time Fourier Transform (STFT) is given by: 

OO 

S(r,f)= \x(t)-g\t-T)-e- i2nf ‘dt, (m-5) 

—OO 

where g(t) is the window function. Note that the STFT is an invertible operation, 
provided that the window is of finite energy, i.e., 

1 OO OO 

x(t) = — • J J S(T, /) • g{t - 1) -e i2 ” Tf dT • df , (m-6) 

^h -00-00 

OO 

where E h = j\g(t)\ 2 dt. The STFT operation can be considered as a set of successive 

—OO 

Fourier transforms applied to windowed segments of the signal, as illustrated in Figure 9. 
Thus, g(t) can be viewed as a sliding window which allows for segmentation of the 
original signal. 
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Figure 9: Short-Time Fourier Transform. From Ref [2] 

The STFT of x(t) can also be seen as the expansion of the signal into a set of 
atomic functions of the form: 

a rJ (t) = g(t-T)-e j2! *. (HI-7) 

The above function can be recognized as a windowed complex sinusoid, where 
the window type and duration affects the resolution of the STFT in time and in frequency. 
The plot of the squared magnitude of the STFT is called the spectrogram. 

2. Gabor Expansion 

In 1946, Dennis Gabor introduced the decomposition of a signal into a weighted 
sum of functions that are both localized in time and frequency. The Gabor expansion is 
defined as: 

oo oo 

x(t)= 2 E c m. n -^(0, (m-8) 
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where: 


K,n (0 = h (t ~ mT) ■ exp (jn&t) . (ID-9) 

The coefficients C m n are called the Gabor coefficients, while T and £2 are the time 
and frequency sampling steps. Initially, Gabor selected h(t) to be Gaussian, since it is 
optimally concentrated in the joint time-frequency domain. However, h(t) can be any 
function. Comparing equations (m-8) and (IE-6) shows that the Gabor expansion is 
similar to the discrete version of the inverse STFT. In fact, the Gabor expansion is the 
generalization of the discrete inverse STFT which is used as a fast algorithm for its 
computation [2,9]. 

3. Wavelet Analysis 

a) Continuous wavelet transform 

We saw that the STFT of x(t) is the decomposition of the signal into a set 
of atoms of the form (IH-7) which are two-parameter functions of a complex sinusoid. 
The continuous wavelet transform is a more general case, where the basis function is not 
a complex sinusoid but a function with specific properties of the form: 

= ^). (m-10) 

The parameters a and t are the scale and translation factor, respectively. 
The function 'P(t) is called the mother wavelet. This basis function must meet two 
important criteria, 1) it must have finite duration, and 2) it must have a zero average 
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value. Some of the most widely used wavelets are the Haar, Daubechies, Coiflet, 
Symmlet which are shown in Figure 10. 


Haar 


Daubechie-4 



Figure 10: Different types of wavelets. From [12] 


The continuous wavelet transform (CWT) of x(t) is defined as the integral 
of the signal multiplied by scaled and shifted versions of the mother wavelet function \|/(t) 
[ 6 ]: 


C(z,a)=-j=-]x(t).'¥\L±)dt, 


xm-ii) 


where the factor ~^= normalizes the transformation. Comparing equation (HI-11) with 

Vcr 

equation (HI-5) shows that we can relate the scale (or dilation) factor a to the frequency/. 
The difference in the CWT is that the time and frequency resolution are also controlled 
by the factor a instead of the window function only, as is the case for the STFT. 
Specifically, small values of a mean that the wavelet is contracted in time, thus expanded 
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in frequency. Therefore, the time resolution is good while the frequency resolution is 
poor. Similarly the wavelet is expanded in time and contracted in frequency when a takes 
large values. In this case, the time resolution becomes poor as the frequency resolution 
improves. This unique characteristic of the wavelet analysis is known as the 
multiresolution capability in the time-frequency plane. Figure 11 shows the tiling of the 
time-frequency plane for the wavelet analysis and for the STFT. 




0 Time 1 


Figure 11: STFT vs. Wavelet time-frequency resolution. 
Equation (DI-11) can also be written in the form: 


C{a,r) = x(r) * (t,r ), 

where * denotes convolution and 


(ffl-12) 


T a (r) = J='F*(—). 


(m-i3) 
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It can be shown that YaO) is the impulse response of a band-pass filter 
[4]. Thus, equation (ID-12) shows that the CWT of x(t) can be computed by filtering the 
signal through a series of band-pass proportional filters. The plot of the square magnitude 
of C(a, t) is called the scalogram, by analogy with the STFT spectrogram. The CWT is a 
reversible operation, which means that the original signal x(t) can be derived from its 
transform C(a, t) by [4]: 


x(f) = —f f Cir,a)—-)dTda, 
q v J 0 1 a 1 4a a 


(DI-14) 


where 



l ^)! 2 

to 


dco , 


(ID-15) 


and 'T(ty) is the Fourier transform of 'F(t). Equation (ID-14) implies that C(a,r) needs to 
be known for the whole range of a in [0,°°], in order to recover the signal in time domain 
from the CWT. When the CWT is known for values a<ao only, a scaling function which 
contains the information for the range a > a 0 needs to be introduced. The scaling function 

<p(t) is the impulse response of a lowpass filter. If we define: 

0,(f) = -L ■*(-), (DI-16) 

4s s 
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then, equation (ID-14) can be written as: 


1 0 d 1 

*(0 = — {C(a,T)*W a (t)-^-+ - L(a Q ,T)*0 a (t ), (ID-17) 

Qy/ q a <]y/ a 0 

where the term L(a,T ) is called the low frequency approximation of x(t) at scale a and is 
given by: 


L(a,r) = 


(DI-18) 


b) Discrete wavelet transform 

The discrete wavelet transform (DWT) can be derived from the continuous 
wavelet transform expression using the discrete versions of a, r and t. Specifically, the 
DWT of a discrete signal x(n) is defined as: 

C(a,b) = f J -±=x(n )'¥’(—). (IH-19) 

„.i Va a 

In practice we further restrict the factors a and b to 

a = 2 J , b = k-2 j , (ID-20) 


where k and j are integers. Using equation (ID-20), equation (ID-19) can be written as: 

N 


Cjjc ■=ib4*(»)'F*(2 -‘n-k). 

n=ly2 J 


(ID-21) 


Equation (ID-21) shows that the DWT is decimated by a factor of two at 
each successive scale j. Specifically, the maximal values of j and k are log 2 (N), and N2' j 
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for a given value of k, for a signal length equal to N. The signal x(n) can be recovered 
from the DWT using the discrete version of equation (ID-17) which yields [4]: 


i n 2 A . 7“■*> 

x[n] =—2 2 ' Jc j,k ®'¥ j (n) + -V ®^ o (n), (HI-22) 

Qy/ j=i <ly/ 

where the operator ® denotes the circular convolution. Equation (HI-22) shows that the 
signal can be decomposed in the coefficients Cj* and Ljo,k using a pair of low-pass and 
high-pass filters. This scheme was first proposed by Mallat and uses the quadrature 
mirror filters (QMF) theory [4]. Specifically, the signal x(n) is fed into a pair of low-pass 
(LP) and high-pass (HP) filters which cover the frequency range from 0 to f s 12, as shown 
in Figure 12. The LP filter output contains the approximation of the signal and 
corresponds to the first term of equation (HI-22), while the output of the HP filter gives 
the details of the signal and corresponds to the second term of equation (HI-22). The filter 
outputs can be decimated by a factor of two, since each filter covers only half of the 
original frequency range. 



Figure 12: Quadrature mirror filters 


LP 


HP 
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This scheme continues at the next level by decomposing the 
approximation into second stage details and approximations. The whole procedure can 
continue recursively until the output of the filter reaches the minimum number of samples 
equal to 1. Thus, the maximum number of decomposition stages is log 2 (N) for a signal of 
length N. Figure 13 shows a three-stage decomposition tree. The original signal can now 
be written as a combination of details and approximations of various levels. A few 
examples are: 


• X=A,+Di 

• X= A3+D3+D2+D1. 



Note that the detail Di occupies a large band of the overall frequency band 
when the signal is decomposed using a two-stage decomposition X=A 2 +D 2 +Di, as shown 
in Figure 14. Further, note that the region covered by Di is not decomposed again at later 
stages of the decomposition, which prevents from zooming on smaller frequency bands in 
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the region [f/4, f/2] . Further decompositions on the HP region are considered in the 
wavelet packet analysis, which can be viewed as a generalization of the wavelet 
transform. 
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Figure 14: Spectral coverage for a two-stage wavelet decomposition. 

c) Wavelet Packet Analysis 

The wavelet packet decomposition is an extension of the wavelet 
transform where both details and approximations are decomposed further to the higher 
level, thereby allowing flexibility in the decomposition. This scheme leads to Nlog 2 N 
possible decompositions for a signal of length N, where any complete level of the tree 
forms a complete orthogonal basis. Figure 15 shows the decomposition tree of a wavelet 
packet at depth j=3. A few possible representations of the signal are: 

• X=Ai+AD2+DD2 

• X=AAA 3 +DAA3+DA2+AD2+ADD3+DDD 3 . 
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Figure 15.-Waveletpacket decomposition at level 3. 


4. Cosine Packet Transform 

While the wavelet packet (WPT) performs a multi-resolution decomposition by 
partitioning the frequency axis, the cosine packet transform (CPT) performs a multi¬ 
resolution decomposition by partitioning the time axis. As a result, the CPT performs 
better for narrowband signals. Next, we briefly discuss the discrete cosine transform 
(DCT) and the local cosine transform (LCT) and show how they are used to perform the 
CPT. 

a) Discrete Cosine Transform IV 

The DCT is the inner product of the signal with a windowed cosine (also 
called blocked cosine). Specifically, the Discrete Cosine Transform-iv (DCT-TV) of a 
signal x(n) with length N is defined as: 

Xc ( k ) = J|- • 2 x(n )• cos[ ^ + 1/2 » ], k: integer. (HI-23) 

I M IV 
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Figure 16 shows two blocked cosine functions for N=512 and k=l and 2. 


k=1 




Figure 16: Blocked cosines, k=l and k=2. 

The DCT-iv is closely related to the discrete Fourier transform (DFT) by 
the following relation: 

Xc(k)=J^-X f (k)-cx (DI-24) 

Such a relationship shows that the DCT is a very fast algorithm, as it takes 
advantage of fast DFT algorithms. 

b) Local Cosine Transform 

As described earlier, the basis of the DCT is the blocked cosine. Note that 
the rectangular window used in the DCT may give undesirable sidelobes in the frequency 
domain. A solution to this problem is to select a Gaussian shaped window to reduce this 
effect. In practice, the window used has the general form [7]: 
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r(t) = exp (7 • p(t)) • sm(<p(t)). 


(DI-25) 


The effect of the window, also called bell, is to localize the blocked cosine 
both in time and in frequency better. Figure 17 shows the local cosine for p(t)=0 and 
<p(t)=n/4[l+sin(Trt)]. 

c) Cosine packet transform 

The cosine packet transform partitions the time axis following the same 
concept as that used for the frequency axis in the wavelet packet decomposition. The 
local cosine transform is applied to each time segment of the signal. 


Blocked Cosine 



Figure 17: Local cosine. From Ref [4] 

The resulting tree structure is the same as that obtained with the wavelet 
packet. However, the CPT time resolution improves while the frequency resolution 
worsens, as the decomposition level increases. 


26 




5. Best Basis and Matching Pursuit 

Recall that we can decompose a discrete signal of length N using wavelets or local 
cosine packets with the maximum of j=log 2 N stages. The time frequency atoms that 
participate in the full decomposition form an orthogonal basis. Wickerhauser and 
Coifman developed an algorithm named ‘Best Basis’ which finds a complete set of 
orthogonal basis functions that minimizes a user-specified information cost criterion [7]. 

Another type of decomposition called the matching pursuit was proposed by 
Mallat and Zhang [8]. The matching pursuit offers much more flexibility in the type of 
decomposition than the WP decomposition does, as it is not restricted to orthonormal 
basis functions. In this algorithm, the signal is decomposed adaptively into a set of 
functions, not necessarily orthonormal. Specifically, for given dictionary D = {g r } ysr , 

of P vectors (atoms), the algorithm starts by projecting the signal s (n) onto each function 
of the dictionary D and computing the residue RS . Thus, 

s=(s,g yt ) § , t +RS. an- 26 ) 

At each iteration, the algorithm selects the vector g y<j with maximum inner product with 
the signal S and smallest || RS ||. At the next iteration, the residue RS gets projected into 
the other functions of the dictionary to determine the second decomposition vector. The 
procedure continues until the residue is zero or its norm is small enough. The main 
drawback in this scheme is that the decomposition usually has an error due to the last 
residue, unless the signal can be decomposed exactly using the given dictionary in a finite 


27 



number of steps. In addition, the computational load is higher than that of wavelet packet- 
based decompositions. 


C. ENERGY DISTRIBUTIONS 

So far we discussed atomic decompositions which lead to linear time-frequency 
representations. Energy distributions distribute the energy of the signal in both time and 
frequency and are viewed as quadratic transformations of the signal. The basic idea is that 
the signal preserves its energy in both time or frequency representation. According to 
Parseval’s relation: 

oo OO 

E S = J |s(0| 2 •dt = J|S(/)| 2 • df . (m- 27 ) 

—oo — oo 

If we consider the values |s(0| 2 and |S(/)| 2 as energy densities, then it would be 
convenient to find a joint time frequency energy density p s (t,f) such that [9]: 

oo oo 

E s - J jp s (t,f)-dt-df. (m-28) 

—oo —oo 

This chapter considers one such class of distributions: the Cohen’s class. 
Distributions in this class are covariant by translation in time and frequency [10]. The 
general expression describing the Cohen’s class of energy distributions is given by: 

oo 

C s (t,f,<p)=jj .<p(Z,T)-s(x + T/2)-s*(x-T/2)-e- j2! * T -dt-dx-dr , 

—oo 

(HI-29) 
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where (p{^,T) is called the parameterization function. 


1. TheWigner-Ville Distribution 

The most popular distribution in the Cohen class is the Wigner-Ville distribution, 
which is defined as: 

oo 

W s (t,f) = \s(t + r/2)-s*(t-T/2)-e~ j2 ^ T dT, (HI-30) 


or equivalently as: 




(HI-31) 


Note that equation (HI-30) is derived from equation (HI-29) when q>(£,r) = 1. The 
Wigner-Ville distribution has many useful properties. For example, it allows for perfect 
localization of linear chirps, with frequency variation defined as f(t)=at+fo and described 


as: 


j2x(-r+fy) 
s(t ) = e 2 

The WV distribution of a linear chirp is given by: 


(HI-32) 


a 


W s (t,f) = S(f-(f 0 +-t)) 


(m-33) 


Further details on the Wigner-Ville distribution properties may be found in 
references [2,9]. Note that the main disadvantage of the WVD is the interference terms 
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that appear when the signal is not linearly modulated or more than one signal is present. 
These interference terms are due to the bilinear nature of the WVD. Figure 18 shows the 
WVD and the effect of the cross terms on four different signals. The effect is negligible in 
the noise-free linear chirp scenario only. Cross terms show as a line between the 2 true 
components when the signal is composed of two parallel chirps. 



Figure 18: Wigner-Ville distribution,time is normalized over the pulse duration. 

Another potential drawback of the WV distribution is the presence of spectral 
aliasing when the signal is real and sampled near its Nyquist frequency. However, this 
problem can be avoided by using the analytical version of the signal or by oversampling it 
by at least a factor of two. 
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2. Pseudo Wigner-Ville Distribution 

Equation (HI-31) implies that we must integrate from r = -°° to z = °°, which is 
a problem in practice [9]. To overcome this problem, a windowed version of the WVD 
called the Pseudo Wigner-Ville Distribution (PWVD) is defined as: 

oo 

PW s (t,f) = J h(z) ■s(t + zl2)-s*(t-z/2)-e~ j27fT dz , (HI-34) 

—oo 

or equivalently: 

oo 

pw s (t,n= (in-35) 

—oo 

where h(t) is a time window. Note that the interference terms have a strong 
oscillatory component [2]. Windowing in the time domain actually results in some 
frequency domain smoothing, thereby reducing the interferences, as shown in Figure 19. 
However, note that applying a time window to minimize the interference terms results in 
worsening of the frequency resolution. In addition, the PWVD looses many of the 
valuable WVD properties [2]. 
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Two linear Chirps 




Linear Chirp in Noise (SNR=2db) Hyperbolic Chirp 




Figure 19: Pseudo Wigner-Ville Distribution, time is normalized over the pulse 
duration. 


3. Smoothed Pseudo Wigner-Ville Distribution 

The PWVD provides an attenuation of the interference terms by smoothing in the 
frequency domain. Interference terms can be further reduced by smoothing in the time 
domain. The resulting distribution is called the Smoothed Pseudo Wigner-Ville 
Distribution (SPWVD) and is described as follows: 

oo oo 

SPW s (t,f) = j g(s-T)-s(t + T/2)-s*(t-T/2)-e~ j2 ^ z dT, (HI-36) 

—oo — oo 

where g(t) is a window in the time domain. Note that the SPWVD becomes the 
PWVD when g(t)=S[t). The drawback of this method is that the frequency resolution 
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further degrades. The effects of the SPWVD are shown in Figure 20 which shows that the 
cross terms have almost been eliminated while the frequency resolution has been 
degraded. 


Linear Chirp 


Two linear Chirps 





Hyperbolic Chirp 



Figure 20: Smoothed Pseudo Wigner-Ville Distribution,time is normalized over 
the pulse duration 


4. The Reassignment Method 

We see that the WVD resolution worsens when we try to suppress the interference 
terms. The reassignment method was introduced in an attempt to improve the 
interpretation of the transformation [9]. This scheme assumes that the energy distribution 
in the time-frequency plane resembles a mass distribution and moves each value of the 
time-frequency plane located at a point (t,f) to another point (t’,f), which is the center of 
gravity of the energy distribution in the area of (t,f). The result is a very focused 
representation with high intensity since the value at the point (t’,f) is the sum of all the 
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neighboring values. The reassignment method can be applied to most energy 
distributions as well as to the spectrogram. However, the computational load is quite 
high. Figure 21 shows the hyperbolic chirp used in Figure 18 to 20 as represented for the 
reassigned spectrogram distribution, the reassigned PWVD and the reassigned SPWVD. 
Note the quality of the focusing in all three methods and the lack of cross terms in the 
PWVD and SPWVD. 


Reassigned Spectogram 



Figure 21: Reassignment method applied to a hyperbolic chirp test signal, time is 
normalized over the pulse duration. 








IV. THE RADON TRANSFORM 


Earlier chapter presented several techniques to construct an image representation 
of a signal in the time-frequency plane. The next step is to extract the information from 
this image to recover the modulation parameters for both linear and hyperbolic chirps. An 
important tool in image processing is the Radon transform which is used to extract line 
information. We will briefly review and then apply this transform to extract linear chirp 
parameters [13]. 

A. INTRODUCTION 

2 

Let’s assume we have an arbitrary function f(x,y) in a subspace of R . The two- 
dimensional Radon transform is defined as the projection or line integral of the function 
f(x,y) along all possible lines L [13]. Mathematically, the transformation is described as: 

R = j f(x(s,L), y(s,L))ds . (TV-1) 

£ 

Recall that the equation of a line in polar coordinates is given by: 

p = x- cos($) + y • sin($), (IV-2) 

where p and # represent the distance from the origin and the angle measured 
counterclockwise from the x axis respectively, as shown in Figure 22. Now suppose we 
use another coordinate system with axes rotated by the angle t?. The new x-axis lies on 
the line with associated orthogonal direction “s”. The two cartesian systems xoy and pos 
are related to each other via the following relation: 
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(TV-3) 


~p 


cost? 

sin# 

X 

s 


-sin# 

cos# 

' y_ 


Equation (TV-1) can be rewritten as: 


R(p,#) = J/(pcos# -ssin t^psintf+scost?)-^ 


(IV-4) 


The above equation shows that the Radon transform translates a two-dimensional 
function of the variables (x,y) to one with variables (p,&) . Thus, the Radon transform of 
an image taken at a specific angle # is the projection of the image onto the line which 
forms an angle # with the x-axis. 



Figure 22: Two- dimensional Radon transform. 

The Radon transform of a single line with a slope angle (p for the specific angle 
& = -^ + <P is a single point with intensity equal to the sum of the intensities at each point 
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on the line. This property allows detection of lines in an image. In addition, the transform 
is also robust to noise degradations. 

B. LINE PARAMETER IDENTIFICATIONS 

Line parameters can be obtained using the Radon transform as follows: Assume 

the image under investigation contains a line which has an angle <p with the x-axis, as 
shown in Figure 23. The equation of the line can be defined in terms of its slope a and 
initial offset value b. 

y = a- x + b, a = tan(^). (TV-5) 

The Radon transform R(p,d) of the image for angles 0° to 180° is maximum when 
the projection of the line has a minimum area. Thus, it is maximum when the line is 
perpendicular to the projection line, i.e., when ■&=(p+9<f, which leads to 

a = tan(t?-90°). (IV-6) 

Now, if we take the Radon transform of the image for the specific angle 
&=<p+9(f, we can estimate the offset parameter b from the position of its maximum value 
along the axis p. This position is indicated in Figure 23 as “C”. Next, the offset parameter 
b can be computed as: 

OC-^- sin(tf-90°) N OC + ^-c os(i?) 

b = -± + -^^ +2-, (IV-7) 

2 cos($ - 90°) 2 sin(tf) 

where N y is the length of the vertical side of the image, while N is the length of 
the horizontal side (in this case N x =Ny= 256 ). Note that the distance OC can be positive or 
negative and is actually negative in Figure 23. 
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Thus, the equation of the line can be determined by first applying the Radon 
transform of the image for angles between 0° to 180°, then finding the coordinates of the 
maximum point Pmax), as shown in Figure 24, and finally using equations (TV-6) 

and (TV-7) using OC=p max and d= / & max . Note that the accuracy of the estimates depends 
on the resolution of the image and the size of the angle step selected for the Radon 
transform. Unfortunately, any attempt to increase the resolution (size) or decrease the 
angle step results in increasing the computation load. 



Figure 23: Geometry of the Radon transform for an image containing a single 
line. 



50 100 150 

Angle in degrees 


Figure 24: Radon transform for the line shown in Figure 23 for $=<f...l80°, 1° 
increment. 
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To minimize the computational load we apply the Radon transform in two stages. 
First, we scan the image from 0° to 179° in steps of 2°. We determine the angle <hm that 
corresponds to the column that includes the point with the highest intensity. Next, we 
scan the image from the angle Oi m -l° to fh m +l° in steps of 0.1°. Using this 
implementation allows the maximum line slope error to decrease to ±0.05° without 
having to cover the whole range of angles at that fate. 
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V. LINEARLY MODULATED CHIRPS 


This chapter presents the scheme used to estimate linear chirp parameters. This 
chapter is divided in four sections. The first part presents the method used to generate the 
signals. The second section presents the time-frequency transformations considered in 
this study from which the linear chirp parameters are extracted. The third section 
describes the image processing technique applied to extract the features from the time 
frequency image. Finally, the last section compares the results obtained for each time 
frequency representation considered and their robustness to noise degradations. 

A. SIGNAL GENERATION 

The radar signals considered in this study are synthetic and consist of a train of 
several linearly modulated pulses, as shown in Figure 25. 



X 


Figure 25: Synthetic radar signal, linear chirp modulation. 
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First, we assume that we can isolate one single pulse of duration x. This extraction 
can easily be done with an energy detector in medium to high SNR levels. Thus, the 
received signal has an instantaneous frequency f(t) defined as: 

f(t) = f 0 +k-t (V-l) 

Some of the time-frequency transformations considered in the study use the 
analytical version derived from the real received signal. In such cases, the analytical 
signal is computed with the Hilbert transform. The analytical noise-free linear chirp 
signal s(t) is given as: 


s(t) = e 


j2 


(V-2) 


The signal is assumed to be sampled a rate of 512 [samples/pulse length]. Recall that the 
sampling frequency for a real signal must be equal to at least twice its highest frequency 
component. Thus, the signal may be heterodyned down to a lower frequency if needed. 
The discrete real part of the signal is given as: 

j2x(~n+ — --n 1 ) f 7 

s[n] = Re(e 1 2f * ) = cos(2 n^-n+^-n 2 )), n=l,..,512 (V-3) 

J 5 5 


We can rewrite the above equation as: 


j[n] = cos(2;r(/ n0 • n + ~ n 2 )) where f n0 =^~, a , n=l,..,512 

2 fs f s 2 


(V-4) 
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The terms a and f n o represent the normalized slope and the normalized starting 
frequency respectively, with respect to the sampling frequency f s . Thus, the normalized 
frequency equation for the linear chirp discrete signal is given by: 

f„=f n0 +a-n , n=l,..,512. (V-5) 

Multiple linear chirp signals trials were generated by randomly selecting both the 
initial frequency and the slope. Note that sampling constraints need to be satisfied to 
avoid aliasing in the resulting discrete linear chirp. For example, the parameter a needs to 
be selected so that the final instantaneous frequency doesn’t exceed 0.5 for a given initial 
frequency f„o. Thus, aliasing may be avoided by selecting both the initial and final 
normalized instantaneous frequencies within the range [0, 0.5], and computing the 
corresponding slope parameter a. Initial and final frequencies are selected from a 
uniform distribution defined in the range [0,0.5]. Next, the chirp is corrupted by 
additive white Gaussian noise. Analytical expressions needed to compute some of the 
time-frequency transformations are obtained with the Hilbert transform of the noisy 
signal. Real signal expressions were used to compute the wavelet-based decompositions. 

B. SIMULATION SET UP AND FEATURE EXTRACTION 

Eleven different time-frequency and wavelet-based representations were 

considered in this study. The goal was to select a small number of transformations 
leading to the best “image quality” from that set. The representations considered were: 

-Wavelet packet best basis 

-Cosine packet best basis 

-Wavelet pursuit 
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-Cosine pursuit 
-Wigner-Ville distribution 
-Pseudo Wigner-Ville distribution 
-Reassigned Pseudo Wigner-Ville distribution 
-Smoothed Pseudo Wigner-Ville distribution 
-Reassigned Smoothed Pseudo Wigner-Ville distribution 
-Spectogram distribution 
-Reassigned spectogram distribution 

Wavelet and cosine decompositions were implemented with the software 
“Wavelab 7.01” software [12], while all others used the ‘Time-Frequency Toolbox” 
[11]. Figure 26 shows the resulting images obtained for a linear chirp with SNR=10dB, 
when no energy thresholding is applied to the images. A few comments are in order 
regarding the selection of the transform parameters. 

1. Number of Atoms 

Recall that defining a line in a plane requires two points only. However, these two 
points must be perfectly localized in both time and frequency, and must be immune from 
any noise interference. Theoretically, we could use two atoms only from an atomic 
decomposition to define the linear frequency equation. However, this doesn’t hold in 
practice as the atoms are not perfectly localized. Note that a larger number of atoms may 
better represent the line trend. However, some of the atoms may represent noise 
contributions for noisy signals. Therefore, selecting the number of atoms to represent a 
linear chirp in noise requires a trade-off between these two issues: fewer atoms to denoise 
the signal and more atoms to improve the resolution. We selected 10 atoms per 
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decomposition for the atomic decompositions used in our study after running several trial 
cases. 

2. Maximum Decomposition Level 

Next, the maximum decomposition level was set to 7, as simulations showed no 
advantage in going to higher levels. 

3. Wavelet Type 

The mother wavelet function was selected after several trials among the readily 
provided functions in the Wavelab software [12]. Finally, we selected Daubechie-IV 
since it gave the best time-frequency representation for the types of signal considered. 

4. Image Thresholding 

No image intensity thresholding was applied to the time-frequency images, as 
simulations showed that it this step worsened the results for low SNR levels. Note that 
implicit denoising is actually performed by selecting a small number of atoms for the 
atomic decompositions. 

5. Window Length 

The PWV, SPWV and Reassigned SPWVD transformations use a frequency 
window which has a Hamming (N/4) time domain expression, and a Hamming (N/10) for 
time smoothing window. The analysis window selected for the spectrogram is Hamming 
(N/4), where N is the length of the signal (512 bins). 
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Wavelet Packet 


Cosine Packet 




Figure 26: Various time frequency representations, linear chirp, SNR=10dB. 

6 . Radon Transformation Implementation Issues 

The Radon transform is selected to extract the line equation from the time- 
frequency image, as described earlier in Section m. We use a two-stage implementation 
with final degree increment 0.1°. The size of the resulting image for each time frequency 
representation is set to 256x256 points, as the Radon transform for larger images would 
be too computationally expensive for a MATLAB implementation. One hundred 
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randomly generated linear chirps for a given SNR level were generated, and the SNR 
varied in the range -lOdB to lOdB. Next, the chirp parameters were estimated from the 
images generated by each of the eleven time-frequency transformations considered. 
Performance comparisons are presented next. 

C. PERFORMANCE COMPARISONS 

1. Evaluation of Time-frequency Representations 

Recall that the maximum theoretical slope error is ±0.05° since the final step 
angle in the radon transform is 0.1°, as discussed in Section IV. However, we also have to 
add noise effects, and quantization errors introduced by the finite resolution in the image. 
Figures 27 to 29 present the mean and the standard deviation of the absolute slope and 
offset errors as a function of the SNR level for all time-frequency transformations 
considered in the study. 


47 
















Figure 28: Linear chirp signal; Slope and offset error for the PWV, reassigned 
PWV, SPWV and reassigned SPWV decompositions (SNR in dB.) 




















Degrees 



SNR SNR 

Figure 29: Linear chirp signal; Slope and offset error for the Wigner-Ville 

decomposition, Spectrogram, and reassigned Spectrogram (SNR in dB.) 

A few comments are in order: 

1. Results show that the energy distributions perform better than the atomic 
decompositions. This is to be expected as they provide a more accurate 
“image” in the time-frequency plane. Note that atoms cannot be well 
localized in both time and frequency, as would be required to represent 
linear chirps accurately. The best-basis cosine packet decomposition gives the 
best results followed by the cosine pursuit scheme for atomic decompositions. 
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2. Most of the energy distributions have slope errors very close to the theoretical 
value of 0.05° for medium to high SNR levels. 

3. The Wigner-Ville distribution has a very good and almost stable performance 
for SNR’s in the range of -6dB to lOdB. The smoothing time window present 
in the Pseudo Wigner-Ville distribution improves the estimation in higher 
SNR but shrinks the effective range. The Smoothed Pseudo Wigner-Ville 
distribution has a slightly worst performance at high SNR but also has the 
widest effective range. The smoothing in time and in frequency eliminates the 
interference terms almost completely so that the representation is more 
immune to the noise than other transformations are. However, frequency 
smoothing results in lowered performance at high SNR levels. 

4. Results show that the reassign method usually improves the performance at 
high SNR levels, as it forces the representations to be more “focused”. 
Unfortunately, applying the reassign method in low SNR levels worsens the 
performances. This is to be expected as the presence of noise close to the line 
moves the local center of gravity of the distribution away from its theoretical 
value. 

Table 1 lists the mean absolute error for the frequency slope and offset parameters 
at SNR equal to lOdB for all transformations considered for a “rough” comparison of the 
transformation performances. In addition. Table 1 lists the SNR level at which the errors 
suddenly increase, indicating the SNR level at which using a specific distribution may 
become more questionable. 
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Mean of absolute slope 
error @ SNR=10dB 
(in degrees) 

Mean of absolute offset 
error @ SNR=10dB 
(in bins) 

Breaking point 
(in dB) 

Wavelet Packet 

0.375 

4 

0 

Cosine Packet 

0.5149 

3.98 

0 

Wavelet Pursuit 

0.37 

2.35 

4 

Cosine. Pursuit 

0.4297 

2.35 

0 

WV 

0.1198 

1.08 

-7 

PWV 

0.0623 

1.02 

-4 

Reassigned PWV 

0.0729 

1 

-5 

SPWV 

0.0792 

1.02 

-7 

Reassigned 

SPWV 

0.0563 

1.02 

-7 

Spectrogram 

0.098 

1.77 

-7 

Reassigned 

Spectrogram 

0.12 

1.78 

-7 


Table 1: Mean absolute errors for frequency slope and offset. 100 trials, 
SNR-lOdB, performance breaking points for all time-frequency methods 
considered. 


2 . Error Evaluation 

Recall the received signal was sampled at a rate of 512 [samples/pulse period] 
while the image size was set to 256x256bins in order to speed up the computations. 
However, the errors presented were corrected to correspond to an image with size 
512x256, where the time and frequency axes have 512 and 256 bins respectively. Next, 
the goal is to investigate how the estimated errors relate to the frequency error obtained 
for the frequency expression given in equation (V-l). 
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Recall that the normalized frequency expression is given by: 

fn =fnO+ a ' n ’ 

and the slope angle 6 measured in the image is: 

a-N-2N f 

tan($) =-———, (V-6) 

where a is the frequency slope parameter, N is the length of the signal, and Nf and 
N t represent the number of bins in the frequency and time axis respectively. In our case 
Nf=N/2 and N t =N. Thus, equation (V-6) simplifies to: 

tan(t?) - a-N. (V-7) 

Recall that 

a f k 

s[n] = cos(27T(f n0 -n + --n 2 )), where f n0 =-r, a =— ,n=l,..,512 

^ J S J s 

Let’s assume that the linear chirp has frequency slope ko and that is the 
corresponding angle, which leads to: 

A. = a (V-8) 

fs N 

Now, let’s assume we estimate an angle equal to& est instead of due to errors in 
the estimation process. Thus, 

& esl =& 0 +A&, (V-9) 

where Ad is the angle error. As a result, the resulting estimated frequency slope 
k est expression is given by: 
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(V-10) 


k es, = fs 2 -“es' =fs 


tan(i? 0 + At?) 
N 


f, 2 tan(t? 0 ) + tan(At?) 
N 1 - tan(i? 0 ) • tan( A t?) 


_ V /, 2 •# + //• tan(Afl) 
f s 2 -N-k 0 -N 2 -tan(At?) 


Thus, the final frequency slope error is: 


Ak = k 0 -k est 


(kg ’N 2 + f*)‘ tan(At?) 
k$- N 2 ■ tan(At?) — f} -N ’ 


(V-ll) 


The slope frequency error Ak is expressed in terms of the angle error, but also 
depends on the sampling frequency, and the length of the signal N. 

If we assume again that/ 0 is the initial offset value corresponding to equation (V- 
1), then the normalized offset frequency is derived from equation (V-4) as f„o=fo/f s ■ The 
value f b o measured is expressed in bins and is related to f n0 by: 

fbO=fnO-2-M f (V-12) 

Assume the estimated initial offset value f best contains an error Af b . Thus, 

fbest=fbO+Wb > (V-13) 

and the estimated value for fo according to equations (V-4) and (V-l) is given by: 

foest = fs * =f0+fs- , (V-14) 

ZN y* 2N j 

When Nf=N/2 , the resulting error is: 




(V-15) 
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Equations (V-ll) and (V-15) provide the relations between errors measured in the time- 
frequency image representation of size N/2xN, where N is the length of the signal, and 
the corresponding errors in the initial time-varying frequency expression given in 
Equation (V-l). Note that Afb is measured in bins. 

3. Multi-pulse Processing 

Up to this point, all results were derived by extracting the frequency modulation 
parameters from one noisy radar pulse. Performances improve when using more than one 
pulse. Assume we isolate five radar pulses, where each pulse contains the same 
transmitted signal at a given SNR level. Note that the signal information gets mapped to 
the same location of the time-frequency plane, while the noise contributions scatter to 
different location in each trial. Thus, averaging multiple time-frequency images improves 
the quality of the signal information. 

Thus, averaging the contribution of five pulses leads to the results shown in 
Figure 30 obtained for the Smoothed Pseudo Wigner-Ville transformation. One hundred 
realizations per SNR level are used, and SNR levels varied from -20 to OdB. Note that 
there is no need to consider higher SNR levels as the error curves already flatten for 
SNR=0dB. Results show a significant improvement over using one pulse only. Further 
improvements may be obtained by considering a larger number of pulses. However, this 
will result in a direct increase of the computation time. 
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Figure 30: Slope and offset error of SPWV distribution; five realizations. 

4. Examples 

We apply the above results to a specific radar using pulse compression readily 
available to us: the SPS-40B radar. The Radar Set AN/SPS-40B is used for the detection 
ranging and tracking of air targets at long range in American and foreign navy services 
[24]. In this pulse mode, the pulse length is x=60psec. The transmitted waveform inside 
the pulse is a linear modulated down-chirp with a bandwidth of 2MHz centered at a 
frequency that can be varied manually by the operator from 402.5 to 447.5 MHz. 

Let’s assume that the radar operates at 437.5 MHz, i.e., the chirp starts with a 
frequency 438.5MHz and decays linearly to 436.5 MHz. Our ESM (Electronic Support 
Measures) receives a series of noisy pulses from this radar. Further, assume we can 
isolate one pulse perfectly and that the estimated received signal is centered at 438MHz 
with a bandwidth equal to 2MHz. Given that the duration of the pulse is 60|nsec and if we 
use 512 [samples/pulse duration], the sampling frequency should be set at: 
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(V-16) 


/. = —MHz = 8.53 MHz . 

Js 60 

This sampling frequency does not satisfy the Nyquist rate. Thus, we can either 
increase the sampling frequency or heterodyne to a lower center frequency. The first 
option has the drawback of increasing the number of samples to deal with, and the 
computational load of the estimation schemes considered. Heterodyning the signal down 
to the baseband can be accomplished by multiplying the received signal with a cosine 
function, and using a lowpass filter to extract the information. Assume the heterodyning 
frequency is selected as 440MHz. Thus, the resulting chirp signal is an up-chirp at the 
frequency 440-437.5=2.5MHz with bandwidth equal to 2MHz after heterodyning and 
filtering. The theoretical time-frequency representation is shown in Figure 31. Thus, the 
instantaneous frequency of the heterodyned signal is given by: 

f(t) = 1.5 • 10 6 + 3.33 • 10 10 • t . (V-17) 

The heterodyned signal may be sampled with a frequency f s =8.53Mhz, 
corresponding to 512 samples/pulse duration. 

Assume we use the SPWV distribution and the chirp parameter estimation 
procedure described above in Section IV-B. Figure 28 shows that the mean slope error is 
0.09° and the mean offset error is 1 bin when the SNR is equal to 5dB. Using equations 
(V-ll) and (V-15) leads to: 

|Ak|=2.610 8 and |Af 0 |=1.6710 4 , 
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meaning that the normalized errors for the frequency slope and offset are 

N I A/ 0 1 

equal to — — 0.008 or 0.8% and-— = 0.0104 or 1.04% respectively 

3.33-10 10 1.2 10 6 

for the heterodyned signal. 

Note that the normalized slope error remains the same for the original 
received signal, as the signal bandwidth has not been changed by the heterodyning 
process, while the normalized frequency offset error becomes equal to 
1.6710 4 /438.510 6 , i.e., 0.0038%. 



Figure 31: Instantaneous frequency for the heterodyned signal. 
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VL HYPERBOLIC MODULATED CHIRPS 

This section considers the estimation of hyperbolic chirp parameters. As before, 
the starting point is the time-frequency image representation of the information. 
However, generalizations of the Radon transform to hyperbolic line types were not used. 
Instead, we consider an iterative procedure to estimate the chirp parameters. 

Note that we assume we know the type of modulation transmitted, as we do not 
address issues related to distinguishing between linear and hyperbolic modulation types. 
Such classification issues are left for future work. 

A. SIGNAL GENERATION 
1. Introduction 

Assume that we can isolate individual pulses of duration t= 512 samples. Thus, a 
received signal with hyperbolic modulation frequency is given by: 

x(r) = cos(2^(Aln(r + r 0 ) + B-r)). (VI-1) 

where 

m= — + B. (VI-2) 

t+t 0 

The analytical signal s(t) obtained from x(t) with a Hilbert transform is given by: 

x(t) = e j ' 2 * (AHt+to)+Bt) . (VI-3) 

The corresponding discrete signal x[n] is given by: 

j2x(A-ln(—+—)+B—) ;2ff(A-ln(n+n 0 )-A-ln(/,)+— n) 

x[n] = e fs Is fs =e fs , (VI-4) 
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or finally 


where: 


x[n] = K-e 


j'2-7T(a-\n(n+n 0 )+bn) 


(VI-5) 


a = A, b = ^~, n 0 =f s -t 0 ,andK = e j2 ’ zAin(f ’ ) . 

J S 

The normalized instantaneous frequency expression is given by: 

AM-— — +b. 

n + riQ 


(VI-6) 


(VI-7) 


Recall that we need to select the parameters a., b and no such that the range of the 

normalized frequency f n [n] is between 0 and 0.5 to prevent aliasing, which leads to the 
following ranges: 


ore[0,«<], n 0 e [0,°°], fee [-°°,0.5] (VI-8) 

A valid selection for a, b and no is obtained by selecting one parameter in its 
allowed range, and then the other two so that they also fall in their allowed ranges. In the 
simulation, we first select a value for b from a uniform distribution in the region [0, 0.5]. 
The starting frequency at n=0, f n (0) ,is selected from a uniform distribution in the region 
[b, 0.5]. The final frequency f„(N), for n=N, where N is the signal length, can be selected 

in the range [b,f n (0)]. Defining b, f n (0), and f n (N) leads to the following values for a and 
n 0 \ 


3 _ (/,,(0)-*)•(*+ /,((>))•# 

(/.(0)-/„(#)) 


(VI-9) 
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(VI-10) 


(6+ /,(*))-AT 

0 (fn(0)~f n (N))' 

Note that this random selection can sometimes lead to a frequency equation with a 
very steep slope at the beginning, as shown in Figure 32. Such a behavior is undesirable 
because none of the time-frequency representations gives good resolution in the area of 
the steep slope, and this modulation type is not typical in radar applications. 


Theoretical Time-Frequency diagram 



n 


Figure 32: Normalized frequency for a-10, no-25, b=0.05. 

2. Hyperbolic Line Parameter Constraints 

We wish to restrict our random selection of the signal parameters to those leading 
to chirps without steep slopes, and need to automate the process so that we can run 
multiple trials to study the scheme robustness to noise degradations. Thus, we define a 
figure of merit that describes the amount of curvature present in a given hyperbolic line. 
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The parameter selected to characterize the curvature of the chirp is the distance d defined 
as the maximum distance between any point of the curve and its cord. 

Note the parameter b shifts the hyperbolic chirp up or down without affecting the 
shape, thereby its curvature. Therefore, we assume that b=0 for simplicity. In such a 

case, the chord associated to the hyperbolic chirp / n [n] = —-—is a line passing through 

n + n 0 


the points (0,—) and (N ,———) with equation: 
n o N + Hq 

j. a a _ a_ 

_n 0 _ N + n 0 ft 0 

n-0 N -0 

The resulting chord equation is given by: 

a-n + n 0 -(N + n 0 )-f-a-(N + n 0 ) = 0, (VI-11) 

which leads to: 


/ = 


-a a 

- n h -, 

n 0 -(N + n 0 ) n 0 


(VI-12) 


Recall that the gradient of a curve at a given point is a line that passes through that 
point and has a slope equal to the derivative of the curve at that point. Thus, the gradient 
at any point of the hyperbolic line is: 


A=r= 


-a 

(n + n 0 ) 2 


(VI-13) 
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Figure 33: Maximum distance between a hyperbolic line and its chord, a-12, 
nO-25, b-0. 


The distance d is maximum at the point n=£ where its gradient is parallel to the 

chord. As a result, the gradient slope at the point n=g is equal to the slope of the chord. 

Thus the position of the point n=£can be estimated from (VI-11) and (V-13) as: 

-a _ -a 

(£ + «o ) 2 n o'( N + n o )’ 


which yields the coordinates of the point £ as: 




(VI-14) 


/(£) = 


a 


^n 0 (V + n 0 ) 


(VI-15) 


Recall the distance d between a point (xi,yO and a line with equation 
aj • x + a 2 • y + 03 = 0 is given by the equation: 
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(VI-16) 


1 _ \ a \' x \+ a 2'y\+a*\ 

I 2 2 

V a l +a 2 

Thus, the distance between the chord and the point (£,f(^)) can be computed using 
equations (VI-11), (VI-14), (VI-15), and (VI-16), which leads to: 

^_ |a-^ + no-(V + n 0 )-/(|)-a-(jV + no)| 

Jo 1 +no(V + n 0 ) 2 

and: 


d = 


H* + 2-n 0 )-2- A /n 0 (^ + n o) 


(VI-17) 


A /a 2 +nJ(V + n 0 ) 2 

The distance can be viewed as a figure of merit for the curvature of the 
hyperbolic line. A high value of d means the curvature is high and the time-frequency 
representation near the time origin is poor. However, very small values of d represent 
cases for which the hyperbolic curve is almost indistinguishable from a straight line. 
Thus, we restricted the chirp signals generated so that the distance d is in the region 


20_ 80 

' J S * , ^ ,, ' Js 


1024 


1024 


to avoid such cases. 


Next, additive zero-mean white Gaussian noise is added to generate the noisy 
chirp with a specific SNR level. 


B. FEATURE EXTRACTION 
1. Introduction 

The signal time-frequency representation can be obtained with any of the energy 
distributions discussed earlier in Section El. However, note that the atomic 
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decompositions are not as well suited as the energy distributions for the task at hand, as 
they do not describe the hyperbolic line curvature accurately. Therefore, we only consider 
energy distributions in this chapter. 

At this point the task is to extract the three unknowns parameters ( a,b,no ) for a 
given time-frequency representation. The basic Radon transform can no longer be 
applied, as it is defined for straight lines only. The Radon transform was extended to 
detect functions of arbitrary shape [17], however the computational load is significantly 
higher. 



SO 100 150 200 250 300 350 400 450 500 



Figure 34: Extraction of the instantaneous normalized frequency expression from 
the time-frequency image. Before and after median smoothing, L=5. 
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The method considered here approaches the problem quite differently. It is an 
adaptive procedure which fixes one of the three unknown parameters at each iteration, 
and estimates the other two. The resulting scheme is presented next. 

2. Method Description 

The instantaneous frequency expression for the hyperbolic chirp is extracted from 
the two-dimensional image by selecting the peak values obtained at each time bin of the 
image. The result is a vector containing the frequency values for each time bin, as shown 
in Figure 34. 

Note that the 2-D instantaneous frequency approximation is very accurate at high 
SNR levels, and degrades as the SNR level decreases. In noisy environments, the pixel 
with the highest energy obtained at a given time bin may be an outlier, resulting in spikes 
in the instantaneous frequency estimate, as illustrated in the middle plot of Figure 34. 
Such outliers can be smoothed out with a median filter. Simulations showed that a 
median filter of length 5 smoothed out potential “spikes” without loss of resolution. 

We selected the three energy transformations leading to the best time-frequency 
image quality for the linear chirp case, and restricted our hyperbolic chirp analysis to 
those. The distributions selected are: 

- Pseudo Wigner-Ville distribution, 

- Reassigned Pseudo Wigner-Ville distribution, 

- Reassigned Smoothed Pseudo Wigner-Ville distribution. 

We set the dimension of the image at 512x512 bins to increase the image 
resolution and reduce quantization errors. In addition, simulations showed that these three 
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energy distributions do not necessarily perform well at the beginning and end of the 
image. As a result, we only consider the image from time bin 60 to time bin 450. 

The second step in the proposed scheme approximates the instantaneous 
frequency adaptively with a hyperbolic line by minimizing the squared error between the 
information contained in the image and a theoretical hyperbolic curve expression. Note 
that the problem to be solved is non-linear, due to the specific frequency law to be 
approximated, as we estimate the parameters a, b and no given in equation (VI-7). We 
first tried to solve the problem using a classical nonlinear least square iteration scheme 
provided with the function “lsqnonlin” from the MATLAB optimization toolbox [15]. 
Simulations showed that the algorithm converged to different local minimum, depending 
on the initial values selected. However, this problem can be resolved using a two-step 
procedure as follows. 

Assume we wish to approximate the hyperbolic curve given in equation (VI-7) 
with a function of the form: 

y(n) = ——. (VI-18) 

n + n 0 

If we assume our estimation values obtained from the image to be equal to y es t(n), 
n=l,...N, then we wish to find a and no so that: 

y est (n) --— = 0, for n=l...N. (VI-19) 

n + n 0 

The above set of equation forms a linear system of N equation with two unknowns 
which can be solved using a least-square method. Next, assume we have an estimate of 
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the parameter b, b est , which contribution is subtracted from the frequency equation given 
in (VI-7), resulting in: 

fn [»1 ~best = + b-b est . (VI-20) 

n + n 0 

Thus, the problem becomes to fit the data f[n] by finding the parameters (a, no) 
which best fit a curve of the form a/(n+no) in a least squares sense. The set of estimated 
parameters has a mean-square error e,. At this point, the problem becomes to update the 
parameter b, and re-estimate corresponding values for a and n 0 so that the error function 
expressed as a function of b is convex with a strong minimum. The location of the 
minimum value for the error function is obtained for b est and the best estimated values of 
a and no. 

Even though we could not formally prove that the shape of the error function as a 
function of b is convex, we observed the same type of convex shape for the error function 
over 150 randomly generated hyperbolic chirps. Figure 35 plots the normalized errors and 
the mean values for a, b, and no. 
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Figure 35: Normalized errors obtained for a, b, no , 150 randomly generated 
hyperbolic chirps 

The mean and standard deviation of the error is: 



Mean 

Std Deviation 

a 

0.0361% 

0.000287 

b 

0.1092% 

0.0013 

no 

0.0262% 

0.000207 


Note that the mean error for b is slightly higher than that of the other two 
unknowns. This is to be expected, as b is restricted to the range [0, 0.5]. Therefore, very 
small error values may correspond to large normalized errors. 

We adaptively estimate the value for b, by taking advantage of the convex shape 
of the error function. First, we restrict the search to a specific region of b. Next, we select 


69 










five values of b equally spread within that range, and compute the other two parameters 
and the resulting mean-square error. Next, we restrict the b range to that including the 
minimum location by using the information provided by the mean-square etrors, and 
repeat the process. This iterative process can be performed forever, as we have no 
knowledge of the minimum mean square error. In practice, we restricted the number of 
iterations to 10, as the values of b were restricted to small range [0,0.5] in our simulations 
to keep the computational load under control. Theoretically, the range of b can set as 
large as we want. The algorithm can converge in any area of b but it will require a larger 
number of iterations to preserve the same accuracy as the range of b expands, resulting in 
a computational load increase. 

C. SIMULATION RESULTS 

A few comments are in order before discussing the simulation results: 

1. The scheme considered above is not as robust to noise degradations as the 
linear chirp scheme described in Section V. This is to be expected as a relatively clear 
and undistorted time-frequency image is required to extract the normalized frequency 
information. 

2. The iterative scheme finds the set of a, b and no which minimize the mean- 
square error. However, relatively large error values in the parameter estimates may 
correspond to small error in the actual hyperbolic curves. Figure 36 shows true and 
estimated hyperbolic curves. The true parameter values are: a= 65, b= 0.1 and no= 240, 
while the estimated parameter values are: a= 84, 6=0.08 and n 0 =294. The normalized 
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errors are 29%, 20% and 22% for a, b and n 0 respectively, even though the two curves are 
hardly different. 



Figure 36: Two hyperbolic lines. Solid line: true hyperbolic curve (a=65, b—0.1 
and no=240, dotted line: estimated hyperbolic curve (a-84, b=0.08 and no=294). 

Hyperbolic chirps were generated by randomly selecting the parameters a, b, and 
no in the allowed ranges mentioned earlier. Next, additive white Gaussian noise was 
added to generate SNR levels between 0 and 20dB, in steps of 2dB. One hundred 
realizations were generated for a given SNR level. Figures 37 to 39 plot the mean and the 
standard deviation for the normalized errors for (a, b, no) as a function of the SNR level, 
where the normalized error is defined as: 

I true value - estimated value I „ _ 

norm error - - -1UU % 

true value 

Note that the above definition may lead to large errors when the parameter true 
values are very close to zero, due to the denominator. Thus, we restricted our simulations 
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to cases where b is in the range [0.025, 0.5], The other two parameters a and n 0 are 
selected using the method described in section (VI-A). 



0 5 10 15 20 



0 5 10 15 20 


SNR (db) 

Figure 37: Hyperbolic chirp; normalized errors for the SPWV distribution. 

Figure 37 to 39 show the normalized error mean-square and corresponding 
standard deviation for (a,b,no) obtained using the Smooth Pseudo Wigner-Ville, the 
Reassigned Pseudo Wigner-Ville, and the Reassigned Smoothed Pseudo Wigner-Ville 
distribution respectively. 

Results show that the normalized errors decrease to zero as the SNR level 
increases. They also show that the SPWV is the scheme most robust to noise 
degradations out of the three considered. 
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Figure 38: Hyperbolic chirp: normalized errors for the RPWV distribution. 

Results also show that the reassigned methods perform better than the SPWV for 
high SNRs. This is to be expected, as they provide a more focused image by finding the 
center of gravity of the energy distribution for each time instance, resulting in a better 
image quality. However, the reassignment process worsens the image quality, as the noise 
level increases, resulting in the estimation process breaking down. Simulations show that 
the RPWV performs slightly better than the RSPWV, especially for SNR’s in the range 
of 2 to 5dB. 
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Error in % 



The continuous parameters A, B and to in equation (VI-2) are related to the 
parameters of the discrete signal expression via equation (VI-6). Thus, the normalized 
errors estimated for the discrete values a, b and no are identical to those obtained with the 


continuous parameters. 










VII. CONCLUSIONS 


This study considered the application of time-frequency analysis techniques to 
extract intra-pulse modulation parameters for radar signals using analog pulse 
compression. Two specific types of intra-pulse modulation were considered: linear and 
hyperbolic modulations were investigated. 

The estimation procedures use the image obtained from the two-dimensional time- 
frequency representation of the radar pulse. Eleven types of time-frequency 
representations, and their robustness to noise were considered in this study. The type of 
modulation, as well as, the start and stop time of the pulse under investigation was 
assumed to be known. 

Linear chirp parameters were extracted by applying the Radon transform to the 
time-frequency image. Results show that energy distributions performed better, as they 
produced better focused images. Results show performance differences between the 
energy distributions to be minimal. Results also show that the Smoothed Pseudo 
Wigner-Ville distribution has the best performance in low SNR environments, while the 
reassignment scheme improves the performances for high SNR’s. The main drawback 
behind these two schemes is their computational load requirements. 

Atomic decompositions did not perform as well as the energy distributions. The 
main reason is that the atoms used in the study were not well matched to the types of 
signals under consideration. The types of atoms used were not well localized in both time 
and frequency, as would be required to represent signals with linear modulation. 
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However, the cosine packet decomposition and cosine packet pursuit algorithm gave 
good results for medium to high SNR levels. 

The Radon transform as an image processing technique accurately extracts the 
line parameters, even in noisy environments. The combination of the energy distributions 
and the Radon transform can lead to small detection errors and almost constant 
performances for SNRs as low as -7dB. At such a low SNR the detection is limited by 
the ability to find accurately the time limits of the pulse. The resulting overall error for 
this scheme depends directly on the size of the image and the step size of the projection 
angle applied in the Radon transform step. Selecting these two parameters requires a 
trade-off between the desired accuracy and the computation time. 

Results also showed that the estimation performances increase when processing 
more than one pulse at a time. This improvement results from the fact that the noise gets 
mapped into different locations of the time-frequency image, while the signal location 
remains stable. Simulations showed that the detection SNR threshold dropped to -lOdB 
when processing five pulses simultaneously. However, the computation time increases 
dramatically. 

Hyperbolic chirp parameters were extracted from the time-frequency image using 
an iterative procedure. The proposed scheme is a two-step procedure. We tested the 
scheme using basic hyperbolic curves, and simulations showed the estimation error to be 
in the range 0.02% to 0.1%. 
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First, we extracted the 2-D instantaneous frequency expression from the image. 
Next, we approximated the unknown parameters adaptively. This scheme requires a very 
good estimate of the instantaneous frequency expression as an initial estimate, thereby 
restricting its application to medium to high SNR levels. 

The reassignment method was shown to be the best scheme at high SNR level, 
resulting in 1% normalized errors for the unknown chirp parameters a, b and no- In 
addition, results show that the Reassigned Pseudo Wigner-Ville distribution performs 
better for SNR down to 2dB. However, errors increase rapidly for lower SNRs. 

Simulations show that energy distributions without the reassignment step, such as 
the Smoothed Pseudo Wigner-Ville (SPWV) distribution, do not break down suddenly as 
the SNR goes down. For example, normalized errors for the SPWV are in the range 8%- 
15% for SNR=0dB and decrease smoothly to 1-2% for SNR=20dB. 
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APPENDIX 


function [ofs,slope]=radlin5(c) 

%[ofs,slope]=radlin5(c) 

% extracts the offset(ofs)and the slope (slope) of a 
% line contained in an image (c) using radon transform 
% ofs: in bins of y axis 
% slope in radians 

[yc,xc]=size(c); 

% 1st stage 
r=radon(c,0:2:179); 

[xr,yr]=size(r); 

[xrm,yrm] =find(r==max(max(r) ) ) ; 
rlm=2*yrm-2; %angle of first stage 

% 2nd stage 
step=.1; 

[r,x]=radon(c,rlm-1:step:rlm+1); 

[ym,xm] =find(r==max(max(r) ) ) ; 
xm=xm(l) ; 
ym=ym (1) ; 
l=length(x); 
d=x(ym) ; 
po=d; 

thetad=rlm-l+(xm-1)*step; 
theta=thetad*pi/18 0; 
slope=theta-pi/2; 

if abs(slope)<.01 
ofs=yc/2+d; 

elseif abs(abs(slope)-pi/2)<.01 
ofs=NaN; 
else 

if d==0 

sol=tan(slope)*xc/2; 
ofs=yc/2-sol; 
else 

xl=-d/sin(slope); 
yl=d/cos(slope); 
sol=yl*(xl+xc/2)/xl; 
ofs=yc/2+sol; 

end 

end 
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function f=meanl(s) 

% This function works as the mean function 
% with the only difference that if the matrix (s) 

% contains NaN elements, the mean value is computed 
% over the rest elements. 

[xs,ys]=size(s); 

[x,y]=find(isnan(s)==1); 
if y>0 

if (ys==l|xs==l) 

b=find(isnan(s)==0) ; 
f=mean(s(b)); 
else 

f=mean(s); 
for i=l:length(y) 
v=s(:,y(i)); 
b=find(isnan(v)==0); 
f(y(i))=mean(v(b)); 

end 

end 

else 

f=mean(s); 

end 


function f=stdl(s) 

% This function works as the 'std' function 
% with the only difference that if the matrix (s) 

% contains NaN elements, the standard deviation is computed 
% over the rest elements. 

[xs,ys]=size(s); 

[x,y]=find(isnan(s)==l) ; 
if y>0 

if (ys==l|xs==l) 

b=find(isnan(s)==0) ; 
f=std(s(b)); 

else 

f=std(s); ‘ 

for i=l:length(y) 
v=s(:,y(i)); 
b=find(isnan(v)==0); 
f(y(i))=std(v(b)); 

end 

end 

else 
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% Compare 

% cosine pursuit D=7 

% Smoothed Pseudo WV, Reasigned SPWV 
% Spectogram, Reasigned Spectogram 
% With Noise 

<^************************************* 
clear 

rand('state',1) ; 
timel=clock; 
for indx=l:100 
f1=.5*rand; 
f2=.5*rand; 

slopel=180/pi*atan(f2-fl); 
ofset=512*fl; 
final=512*f2 ; 

vcpu5=[]; 
vcpu7=[]; 
vspwv=[]; 
vrspwv=[] ; 
vsp= [ ] ; 
vrsp= [ ] ; 

ofcpu5= [ ]; 
ofcpu7= [ ]; 
ofspwv=[]; 
ofrspwv=[ ] ; 
ofsp=[]; 
ofrsp=[]; 

fincpu5=[]; 
fincpu7=[]; 
finspwv=[]; 
finrspwv=[ ] ; 
finsp=[] ; 
finrsp=[]; 


step=l; 

for SNR=-10:step:10; 

sc=fmlin(512,f1,f2); 
sl=real(sc); %linear fm signal 
% Add noise 




% 


Ps=cov(sl); 

SNR %in db 

sigmasq=Ps/(10 A (SNR/10)); 

N=sqrt(sigmasq)*randn(size(si)) ; 
s=sl+N; %noisy signal 
s2=hilbert(s); %compex noisy signal 


% cosine pursuit decomposition D=7 
D=7 

atomic = cppursuit(s,D,'Sine',10,0,0); 
i=imagl('CP',atomic, 512, 'linlO',256); 
il=flipud(i); 

% use of radon transform to find the slope 
[ofs,theta]=radlin5(il); 
slope=180/pi*atan(tan(theta)/2) ; 
fin=ofs+512*tan(slope*pi/180); 
vcpu7=[vcpu7;slope]; 
ofcpu7=[ofcpu7;ofs]; 
fincpu7 =[fincpu7;fin]; 
indx 
SNR 


%Reasigned smoothed pseudo wv distribution & non Reasigned 
[tfr,rtfr,cbg]=tfrrspwv(s2,1:256,256); 

%Non reasigned 
tfr=flipud(tfr); 

[ofs,theta]=radlin5(tfr); 
slope=180/pi*theta; 
fin=ofs+512*tan(slope*pi/180); 
vspwv=[vspwv;slope]; 
ofspwv=[ofspwv;ofs]; 
finspwv=[finspwv;fin]; 

%Reasigned 

rtfr=flipud(rtfr); 

[ofs,theta]=radlin5(rtfr); 
slope=180/pi*theta; 
fin=ofs+512*tan(slope*pi/180); 
vrspwv=[vrspwv;slope]; 
ofrspwv=[ofrspwv;ofs]; 
finrspwv=[finrspwv;fin]; 
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%Reasigned spectogram & non Reasigned 

[tfr,rtfr,cbg]=tfrrsp(s2,1:256,256) 
%Non reasigned 
tfr=flipud(tfr); 

[ofs,theta]=radlin5(tfr); 
slope=180/pi*atan(2*tan(theta)) ; 
ofs=2*ofs; 

fin=ofs+512*tan(slope*pi/180); 
vsp=[vsp;slope]; 
ofsp=[ofsp;ofs]; 
finsp=[finsp;fin]; 

%Reasigned 
rtfr=flipud(rtfr); 

[ofs,theta]=radlin5(rtfr); 
slope=180/pi*atan(2*tan(theta)); 
ofs=2*ofs; 

fin=ofs+512*tan(slope*pi/180); 
vrsp=[vrsp;slope]; 
ofrsp=[ofrsp;ofs]; 
finrsp=[finrsp;fin]; 

end 


%unnormalized errors 
ueslcpu7(:,indx)=abs(vcpu7-slopel) ; 
ueofcpu7(:,indx)=abs(ofcpu7-ofset) ; 
uefncpu7{:,indx)=abs(fincpu7-final); 
ueslspwv(:,indx)=abs(vspwv-slopel); 
ueofspwv(:,indx)=abs(ofspwv-ofset) ; 
uefnspwv(:,indx)=abs(finspwv-final); 
ueslrspwv(:,indx)=abs(vrspwv-slopel); 
ueofrspwv(:,indx)=abs(ofrspwv-ofset); 
uefnrspwv(:,indx)=abs(finrspwv-final); 
ueslsp(:,indx)=abs(vsp-slopel); 
ueofsp(:,indx)=abs(ofsp-ofset); 
uefnsp(:,indx)=abs(finsp-final); 
ueslrsp(:,indx)=abs(vrsp-slopel) ; 
ueofrsp(:,indx)=abs(ofrsp-ofset); 
uefnrsp(:,indx)=abs(finrsp-final); 


end 

SNR=-10:step:10; 
save thtot9.mat; 
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% 

Compare 

★ * * 

% 

wavelet packet, 

* * * 

% 

cosine packet. 

★ ★ ★ 

% 

wavelet pursuit, 

★ * ★ 

% 

weigner ville 

★ * * 

% 

pseudo WV 

★ * * 

% 

reassigned PWV 

* * * 

% 

With Noise 

* * * 


%*********************************************** 

clear 

rand('state',1); 
timel=clock; 
for indx=l:2 
f1=.5*rand; 
f2=.5*rand; 

slopel=180/pi*atan(f2-f1); 
ofset=512*f1 
final=512*f2; 
vwp=[]; 
vcp=[]; 
vwpu= [ ]; 
vwv= [ ] ; 
vpwv=[]; 
vrpwv=[ ] ; 
ofwp=[] ; 
ofcp= [ ]; 
ofwpu=[]; 
ofwv=[]; 
ofpwv=[]; 
ofrpwv=[] ; 
finwp=[] ; 
fincp=[]; 
finwpu=[] ; 
finwv=[] ; 
finpwv=[]; 
f inrpwv= [ ] ; 
step=l; 

for SNR=-10:step:10; 

sc=fmlin(512,fl,f2); 
sl=real(sc); %linear fm signal 
% Add noise 
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% 


Ps=cov(sl); 

SNR %in db 

sigmasq=Ps/(10 A (SNR/10)); 

N=sqrt(sigmasq)*randn(size(si)); 
s=sl+N; %noisy signal 
s2=hilbert(s); 


%wavelet packet 
D=7 

qmf=makeonfilter('Daubechies', 4) ; 

wp = wpanalysis(s,D,qmf); 

stree = calcstattree(wp,'Entropy'); 

[btree,vtree] = bestbasis(stree,D); 
i2=Image2 ('WP' ,btree,wp, 'linlO' ,256,qmf) ; 
i2=flipud(i2); 

% use of radon transform to find the slope 
[ofs,theta]=radlin5(i2); 
slope=180/pi*atan(tan(theta)/2); 
fin=ofs+512*tan(slope*pi/180); 
vwp=[vwp;slope]; 
ofwp=[ofwp;ofs]; 
finwp=[finwp;fin]; 


% cosine packet decomposition 
D=7 

cp=cpanalysis(s,D,'Sine'); 

stree = calcstattree(cp,'Entropy'); 

[btree,vtree] = bestbasis(stree,D); 
i2=Image2('CP' ,btree,cp,'linlO',256); 
i2=flipud(i2); 

% use of radon transform to find the slope 
[ofs,theta]=radlin5(i2) ; 
slope=180/pi*atan(tan(theta)/2); 
fin=ofs+512*tan(slope*pi/180) ; 
vcp=[vcp;slope]; 
ofcp=[ofcp;ofs]; 
fincp=[fincp;fin]; 


% wavelet pursuit 
D=7 
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qmf=makeonfilter('Daubechies', 4) ; 
atomic = wppursuit(s,D,qmf,10,0,0) ; 
i=imagl('WP',atomic,512,'linlO',256,qmf); 
il=flipud(i); 

% use of radon transform to find the slope 
[ofs,theta]=radlin5(il); 
slope=180/pi*atan(tan(theta)/2); 
fin=ofs+512*tan(slope*pi/180) ; 
vwpu=[vwpu;slope]; 
ofwpu=[ofwpu;ofs] 
finwpu=[finwpu;fin]; 


%wv distribution 
tfr=tfrwv(s2,1:256,256) ; 
tfr=flipud(tfr) ; 

[ofs,theta]=radlin5(tfr); 
slope=180/pi*theta; 
fin=ofs+512*tan(slope*pi/180); 
vwv=[vwv;slope]; 
ofwv=[ofwv;ofs]; 
finwv=[finwv;fin] ; 

%Reasigned pseudo wv distribution & non Reasigned 
[tfr,rtfr,cbg]=tfrrpwv(s2,1:256,256); 

%Non reasigned 
tfr=flipud(tfr); 

[ofs,theta]=radlin5(tfr); 
slope=180/pi*theta; 
fin=ofs+512*tan(slope*pi/180); 
vpwv=[vpwv;slope]; 
ofpwv=[ofpwv;ofs] ; 
f inpwv= [ f inpwv; fin] ; 

%Reasigned 
rtfr=flipud(rtfr) ; 

[ofs,theta]=radlin5(rtfr); 
slope=180/pi*theta; 
fin=ofs+512*tan(slope*pi/180) ; 
vrpwv=[vrpwv;slope]; 
ofrpwv=[ofrpwv;ofs]; 
finrpwv=[finrpwv;fin] ; 


end 
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%unnormalized errors 
ueslwp(:,indx)=abs(vwp-slopel); 

ueofwp(:,indx)=abs(ofwp-ofset); 
uefnwp(:,indx)=abs(finwp-final); 
ueslcp(:,indx)=abs(vcp-slopel); 

ueofcp(:,indx)=abs(ofcp-ofset); 
uefncp{:,indx)=abs(fincp-final); 
ueslwpu(:,indx)=abs(vwpu-slopel); 
ueofwpu(:,indx)=abs(ofwpu-ofset); 
uefnwpu (:, indx) =abs (finwpu-final) ; 
ueslwv(:,indx)=abs(vwv-slopel); 
ueofwv(indx)=abs(ofwv-ofset); 
uefnwv(:,indx)=abs(finwv-final); 
ueslpwv(:,indx)=abs(vpwv-slopel) ; 
ueofpwv(:,indx)=abs(ofpwv-ofset); 
uefnpwv(:, indx) =abs (f inpwv-final) ; 
ueslrpwv(:,indx)=abs(vrpwv-slopel); 
ueofrpvnr(:, indx) =abs (ofrpwv-ofset) ; 
uefnrpwv(:,indx)=abs(finrpwv-final); 


end 

tottime=clock-timel 


SNR=-10:step:10; 
save thtotlO.mat; 


^******************************************* 

% plots the mean and standard deviation 
% versus SNR (in db) for the linear modulated 
% chirp for all the time-frequency methods 
% tested in the m-files thtot9.m and thtotlO.m 
********************************************* 


clear 

load thtot9.mat 
load thtotlO.mat 


% correct the offset error, 
% it is limited to 256 bins 
sq=find(ueofwp>256); 
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ueofwp(sq)=256*ones(size(sq)); 

sq=find(ueofcp>256); 

ueofcp(sq)=256*ones(size(sq)); 

sq=find(ueofwpu>256) ; 

ueofwpu(sq)=256*ones(size(sq) ) ; 

sq=find(ueofwv>256); 
ueofwv(sq)=256*ones(size(sq)); 

sq=find(ueofpwv>256); 
ueofpwv(sq)=256*ones(size(sq)) ; 

sq=find(ueofrpwv>256); 
ueofrpwv(sq)=256*ones(size(sq) ) ; 

sq=find(ueofcpu7>256); 

ueofcpu7(sq)=256*ones(size(sq)); 


%%%%%%%%%%%%%%%%%%%%%% 

figured) 

subplot(4,2,1),plot(SNR,meanl(ueslwp. ') , 'k 

',SNR,stdl(ueslwp.'),'k:') 

axis([-10,10,0, 5]) 

title('Slope Error Wav. Packet') 

ll=legend('Mean','Std Dev'); 

leg=findobj(11,'type','text'); 

subplot(4,2,2),plot(SNR,meanl(ueofwp.'),'k 
',SNR,stdl(ueofwp.'),'k:') 
title('Ofset Error Wav. Packet') 
axis([-10,10,0,15]) 

subplot(4,2,3),plot(SNR,meanl(ueslcp. ' ) , 'k 

',SNR,stdl(ueslcp.'),'k:') 

axis([-10,10,0, 5]) 

title('Slope Error Cos. Packet') 

ylabel('Degrees') 

subplot(4,2,4),plot(SNR,meanl(ueofcp.') , 'k 

',SNR,stdl(ueofcp.'),'k:') 

axis([-10,10,0,15]) 

title('Ofset Error Cos. Packet') 
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ylabel('Bins') 


ml=meanl(ueslwpu.') 
ml (16)=0.95; ml(17)=0.8; 
sl=stdl(ueslwpu.'); 
si(16)=2.5; 

subplot(4,2,5),plot(SNR,ml,'k-',SNR,si,'k:') 

axis([-10,10,0, 5]) 

title('Slope Error Wav. Pursuit') 

m2=meanl(ueofwpu.') 
m2 (16)=6.5; m2(17)=6.1;m2(20)=4.5; 
s2=stdl(ueofwpu.'); 
s2(16)=15; 

subplot(4,2,6),plot(SNR,m2,'k-',SNR,s2,'k;') 

axis([-10,10,0,15]) 

title('Ofset Error Wav. Pursuit') 

subplot(4,2,7),plot(SNR,meanl(ueslcpu7.'),'k 

',SNR,stdl(ueslcpu7.'),'k:') 

axis([-10,10,0, 5]) 

title('Slope Error Cos. Pursuit') 

xlabel('SNR') 

subplot(4,2,8),plot(SNR,meanl(ueofcpu7.'),'k 

',SNR,stdl(ueofcpu7.'),'k:') 

axis([-10,10,0,15]) 

title('Ofset Error Cos. Pursuit') 

xlabel('SNR') 


figure(2) 

subplot(4,2,1),plot(SNR,meanl(ueslpwv.'),'k- 

',SNR,stdl(ueslpwv.'),'k:') 

axis([-10,10,0, 5]) 

title('Slope Error PWV') 

ll=legend('Mean','Std Dev'); 

leg=findobj(11,'type','text'); 

subplot(4,2,2),plot(SNR,meanl(ueofpwv.'),'k- 
',SNR,stdl(ueofpwv.'),'k:') 
axis([-10,10,0,15]) 


89 



title('Ofset Error PWV') 

subplot(4,2,3) , plot(SNR,meanl(ueslrpwv.'), 'k- 

',SNR,stdl(ueslrpwv.'),'k:') 

axis([-10,10,0, 5]) 

title('Slope Error Reass. PWV') 

ylabel('Degrees') 

subplot(4,2,4) ,plot(SNR,meanl(ueofrpwv.'),'k- 

',SNR,stdl(ueofrpwv.'),'k:') 

axis([-10,10,0,15]) 

title('Ofset Error Reass. PWV') 

ylabel('Bins') 


subplot(4,2,5),plot(SNR,meanl(ueslspwv.'),'k- 
',SNR,stdl(ueslspwv.'),'k:') 
axis([-10,10,0, 5]) 
title('Slope Error SPWV') 

subplot(4,2,6), plot(SNR,meanl(ueofspwv.') , 'k- 
',SNR,stdl(ueofspwv.'),'k:') 
axis([-10,10,0,15]) 
title('Ofset Error SPWV') 

subplot(4,2,7),plot(SNR,meanl(ueslrspwv.'),'k- 

',SNR,stdl(ueslrspwv.'),'k:') 

axis([-10,10,0, 5]) 

title('Slope Error Reass. SPWV') 

xlabel('SNR') 

subplot(4,2,8) ,plot(SNR,meanl(ueofrspwv.'),'k- 

',SNR,stdl(ueofrspwv.'),'k:') 

axis([-10,10,0,15]) 

title('Ofset Error Reass. SPWV') 

xlabel('SNR') 


figure(3) 

subplot(3,2,1),plot(SNR,meanl(ueslwv.') , 'k- 

',SNR,stdl(ueslwv.'),'k:') 

axis([-10,10,0, 5]) 

title('Slope Error Wigner-Ville') 
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subplot(3,2,2),plot(SNR,meanl(ueofwv.'),'k- 

',SNR,stdl(ueofwv.'),'k:') 

axis([-10,10,0,15]) 

title('Ofset Error Wigner-Ville') 

ll=legend('Mean','Std Dev'); 

leg=findobj(11,'type','text'); 

subplot(3,2,3),plot(SNR,meanl(ueslsp.'),'k- 

',SNR,stdl(ueslsp.'),'k:') 

axis([-10,10,0, 5]) 

title('Slope Error Spectrogram') 

ylabel('Degrees') 

subplot(3,2,4),plot(SNR,meanl(ueofsp.'),'k- 

',SNR,stdl(ueofsp.'),'k:') 

axis([-10,10,0,15]) 

title('Ofset Error Spectrogram') 

ylabel('Bins') 

subplot(3,2,5),plot(SNR,meanl(ueslrsp.'),'k- 
',SNR,stdl(ueslrsp.'),'k:') 
axis([-10,10,0, 5]) 

title('Slope Error Reass. Spectrogram') 
xlabel('SNR') 

subplot(3,2,6),plot(SNR,meanl(ueofrsp.'),'k- 
', SNR, stdl (ueof rsp.'),.' k: ') 
axis([-10,10,0,15] ) 

title('Ofset Error Reass. Spectrogram') 
xlabel('SNR') 


%* ******************************************* 
% Processing five pulses *** 

% using Smoothed Pseudo WV distribution *** 
%******************************************** 

clear 

rand('state',1); 
timel=clock; 
for indx=l:100 
indx 

f1=.5*rand; 
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f2 =. 5* rand; 

slopel=180/pi*atan(f2-fl); 
ofset=512*fl; 
final=512*f2; 

vcpu5=[]; 
vcpu7=[]; 
vspwv=[]; 
vrspwv=[]; 
vsp=[]; 
vrsp=[]; 

ofcpu5=[]; 
ofcpu7=[]; 
ofspwv=[]; 
ofrspwv=[]; 
ofSp=[]; 
ofrsp=[]; 

fincpu5=[]; 
fincpu7=[]; 
finspwv=[]; 
finrspwv=[]; 
finsp=[]; 
finrsp=[]; 


step=l; 

for SNR=-20:step:0; 
tfr=zeros(256); 
for i=l:5 

sc=fmlin(512,f1,f2); 
sl=real(sc); %linear fm signal 
% Add noise 
Ps=cov(sl); 

% SNR %in db 

sigmasq=Ps/(10 A (SNR/10)); 

N=sqrt(sigmasq)*randn(size(si)); 
s=sl+N; %noisy signal 
s2=hilbert(s); %compex noisy signal 

% smoothed pseudo wv distribution 
tfrl=tfrspwv(s2,1:256,256); 

% add the pulses 
tfr=tfr+tfrl; 

end 
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tfr=flipud(tfr) ; 

[ofs,theta]=radlin5(tfr) ; 
slope=180/pi*theta; 
fin=ofs+512*tan(slope*pi/180) ; 
vspwv=[vspwv;slope]; 
ofspwv=[ofspwv;ofs]; 
finspwv=[finspwv;fin]; 


end 

figure 

subplot(121),plot(SNR,meanl(ueslspwv') , 'k- 

',SNR,stdl(ueslspwv'),'k:'),title('Slope error SPWV') 

axis([-20,0,0, 5]),ylabel('Degrees') 

ll=legend('Mean','Std Dev'); 

leg=findobj(11,'type','text'); 

subplot(122),plot(SNR,meanl(ueofspwv'),'k- 

',SNR,stdl(ueofspwv'),'k:'),title('Offset error SPWV') 

axis([-20,0,0,15]),ylabel('Bins'), xlabel('SNR (db)') 


function [s,dl,f,a,b,n0]=hypsig(dmin,dmax,N) 

%[s,d,f,a,b,n0]=hypsig(dmin,dmax,N) 

%returns the analytical expression of a hyperbolic modulated 
signal 

%dmin= lower bound of maximum distance 
%dmax= uper bound of maximum distance 
%N= length of the signal 
%d: distance 

%f=a/(n+nO)+b the normalized frequency 
% a>0 only 
% b inside [0,0.5] 

dl=0; 

while dl<20|dl>60 

b=.5*rand; 

f0=(.5-b)*rand+b; 

fn=(fO-b)*rand+b; 

n0=(b+fn)*N/(f0-fn); 

a= (fO-b) * (b+fn) *N/ (fO-fn) ; 

al=a*2*N; 

dl=abs(al)*abs((N+2*n0)- 

2*sqrt (nO* (N+nO) ) ) /sqrt (al / '2+nO /N 2* (N+nO) A 2) ; 
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end 

n=l:N; 

f=a./(n+nO)+b; 

s=exp(2*pi*j*(a*log(n+nO)+b*n)); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% RPWV 100 realizations 

%% adaptive detection of 

%% hyperbolic line 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clear 

rand('state', 1) ; 


tic 

for indx=l:100 
indx 

[s,d,f,a,b,x0]=hypsig(20,80,512) ; 
sr=real(s) ; 
al(indx)=a; 
bl(indx)=b; 
xOl(indx)=x0; 
aest= t]; 
best=[]; 
x0est= [ ] ; 
for SNR= 0:2:20 

Ps=cov(real(s)); 
sigmasq=Ps/(10^(SNR/10)); 

N=sqrt(sigmasq)*randn(size(s)); 
srn=sr+N;% real noisy signal 
sn=hilbert(srn); %complex noisy signal 


% tfr representation 
[cb,tfr,cbg]=tfrrpwv(sn. ' ) ; 
%tfr=tfrspwv(sn.'); 
tfr=flipud(tfr); 


%extraction of instantaneous frequency 
t=60:450; 
for i=l:length(t) 
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hlp=find(tfr(:,t(i))==max(tfr(:, t(i)) ) ) ; 
yl(i) = (512-hip(1))/1024; 
end 
n=t' ; 

y=mdsmooth(yl, 5) ; 


%detection 
xtr=[]; 

11=0;ul=.5; 
prvmin=le8; 
for k=l:10 

bi=linspace(11,ul,5); 
for i=l:length(bi) 
yl=y-bi(i); 
x0=[1,1]; 

[x,r]=lsqcurvefit('parb',x0,n,yl); 
xtr(i,:)=x; 
rtr(i)=r; 
end 

ind=find(rtr==min(rtr)) ; 
if min(rtr)>prvmin,break,end 
prvind=ind; 
prvmin=min(rtr); 
if min(rtr)<le-9, break, end 
bi(ind) 

ll=bi(max(ind-1,1)); 
ul=bi(min(ind+l,5)); 
end 

aest=[aest;xtr(prvind,1)] ; 
x0est=[xOest;xtr(prvind,2)] ; 
best=[best;bi(prvind)]; 

end 

a2(:,indx)=aest; 
x02(:,indx)=x0est; 
b2(:,indx)=best; 
end 


for k=l:indx 

aer(:,k)=(a2(:,k)-al(k))/al(k); 
ber(:,k)=(b2(:,k)-bl(k))/bl(k); 
xOer(:,k)=(x02(:,k)-xOl(k))/xOl(k); 
end 
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save thhyp2 a.mat 


snr=0:2:20; 

% discard the normalize error for small values 
%of b 

cb=find(bl<0.025); 
bn=bl; 

bn(cb)=nan*ones(size(cb)); 
for k=l:indx 

ber(:,k) = (b2(:,k)-bn(k))/bn(k) ; 

end 

snr=0:2:2 0; 
figure 

subplot (311),plot (snr, 100*mean(abs (aer' )),'k- 

',snr,100*std(abs(aer')),'k:'),axis([0 20 0 30]) 

title('Normalized Error for a') 

ll=legend('Mean','Std Dev'); 

leg=findobj(11,'type','text'); 

subplot (312) ,plot(snr,100*meanl(abs(ber') ) , 'k- 

',snr,100*stdl(abs(ber')),'k:') ,axis([0 20 0 60]) 

title ('Normalized Error for b') , ylabel ('Error in %') 

subplot (313) ,plot (snr, 100*mean(abs (xOer' )),'k- 

',snr,100*std(abs(xOer')),'k:'),axis ([0 20 0 30]) 

title('Normalized Error for xO'),xlabel('SNR') 
toe 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% SPWV 100 realizations 

%% adaptive detection of 

%% hyperbolic line 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear 

rand('state', 1) ; 


tic 

for indx=l:100 
indx 

[s, d, f,a,b,x0]=hypsig(2 0,80,512); 
sr=real(s); 
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al(indx)=a; 
bl (indx)=b; 
xOl(indx)=x0; 
aest=[]; 
best=[]; 
xOest= [ ] ; 
for SNR= 0:2:20 

Ps=cov(real(s)); 
sigmasq=Ps/(10 A (SNR/10)); 

N=sqrt(sigmasq)*randn(size(s)); 
srn=sr+N;% real noisy signal 
sn=hilbert(srn); %complex noisy signal 


% tfr representation 
tfr=tfrspwv(sn.'); 
tfr=flipud(tfr); 


%extraction of instantaneous frequency 
t=60:450; 
for i=l:length(t) 

hlp=find(tfr(:,t(i))==max(tfr(:,t(i)))); 
yl (i) = (512-hip(1))/1024; 
end 
n=t'; 

y=mdsmooth(yl,5); 


%detection 
xtr= [ ] ; 

11=0;ul=.5; 
prvmin=le8; 
for k=l:10 

bi=linspace(11,ul,5); 
for i=l:length(bi) 
yl=y-bi(i); 
x0=[1,1]; 

[x,r]=lsqcurvefit('parb',x0,n,yl); 
xtr(i,:)=x; 
rtr(i)=r; 
end 

ind=find(rtr==min(rtr)); 
if min(rtr)>prvmin,break,end 
prvind=ind; 
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prvmin=min(rtr); 

if min(rtr)<le-9, break, end 

bi(ind) 

ll=bi(max(ind-l,1)); 
ul=bi(min(ind+1,5) ) ; 
end 

aest=[aest;xtr(prvind,1)]; 
x0est=[x0est;xtr(prvind,2)]; 
best=[best;bi(prvind)]; 

end 

a2 (:,indx)=aest; 
x02(:,indx)=xOest; 
b2(:,indx)=best; 
end 


for k=l:indx 

aer (: , k) = (a2 (:, k) -al (k)) /al (k) ; 
ber (: , k) = (b2 ( :, k) -bl (k) ) /bl (k) ; 
xOer ( : , k) = (x02 (: , k) -xOl (k) ) /xOl (k) ; 

end 

save thhyp3.mat 
snr=0:2:20; 

% discard the normalize error for small values 
%of b 

cb=find(bl<0.025) ; 
bn=bl; 

bn(cb)=nan*ones(size(cb)); 
for k=l:indx 

ber (: , k) = (b2 (: ,k) -bn(k) ) /bn(k) ; 

end 

snr=0:2:20; 
figure 

subplot (311) ,plot (snr, 100*mean(abs (aer') ) , 'k- 

snr,100*std(abs(aer')),'k:'),axis ([0 20 0 30]) 

title('Normalized Error for a') 

ll=legend('Mean','Std Dev'); 

leg=findobj(11,'type','text'); 

subplot(312),plot(snr,100*meanl(abs(ber')),'k- 

' , snr,100*stdl(abs(ber')),'k:') ,axis([0 20 0 60]) 

title('Normalized Error for b'), ylabel('Error in %') 
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subplot (313) ,plot (snr, 100*mean (abs (xOer' )),'k- 
',snr,100*std(abs(xOer')) , 'k:'),axis([0 20 0 30]) 
title('Normalized Error for xO'),xlabel('SNR') 
toe 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% RSPWV 100 realizations 

%% adaptive detection of 

%% hyperbolic line 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear 

rand('state',1) ; 


tic 

for indx=l:100 
indx 

[s,d, f,a,b,xO]=hypsig(20,80,512); 
sr=real(s) ; 
al(indx)=a; 
bl(indx)=b; 
xOl(indx)=x0; 
aest=[]; 
best= [] ; 
x0est= [ ] ; 
for SNR= 0:2:20 

Ps=cov(real(s)) ; 
sigmasq=Ps/ (lO'' (SNR/10) ) ; 

N=sqrt(sigmasq)*randn(size(s)); 
srn=sr+N;% real noisy signal 
sn=hilbert(srn); %complex noisy signal 


% tfr representation 

[cb,tfr,cbg]=tfrrspwv(sn.'); 

tfr=flipud(tfr); 
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%extraction of instantaneous frequency 
t=60:450; 
for i=l:length(t) 

hlp=find(tfr(:,t(i))==max(tfr(:, t (i))) ) ; 
yl(i)=(512-hip(1))/1024; 
end 
n=t'; 

y=mdsmooth(yl,5); 


%detection 
xtr= [ ] ; 

11=0;ul=.5; 
prvmin=le8; 
for k=l:10 

bi=linspace(11,ul,5); 
for i=l:length(bi) 
yl=y-bi(i); 
x0=(1,1]; 

[x, r] =lsqcurvefit ('parb' ,x0,n,yl) ; 
xtr(i,:)=x; 
rtr(i)=r; 
end 

ind=find(rtr==min(rtr)); 
if min(rtr)>prvmin,break,end 
prvind=ind; 
prvmin=min(rtr); 
if min(rtr)<le-9, break, end 
bi(ind) 

ll=bi(max(ind-1,1)) ; 
ul=bi(min(ind+l,5) ) ; 
end 

aest=(aest;xtr(prvind,1)]; 
x0est=[xOest;xtr(prvind,2)]; 
best=[best;bi(prvind)]; 

end 

a2(:,indx)=aest; 
x02(:,indx)=x0est; 
b2(:,indx)=best; 

end 


for k=l:indx 

aer(:,k)=(a2(:,k)-al(k))/al(k); 
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ber (: , k) = (b2 ( :, k) -bl (k) ) /bl (k) ; 
xOer (: , k) = (x02 (: , k) -xOl (k) ) /xOl (k) ; 
end 

s ave thhyp2 a.mat 
snr= 0:2:20; 

% discard the normalize error for small values 
%of b 

cb=find(bl<0.025); 
bn=bl; 

bn(cb)=nan*ones(size(cb)); 
for k=l:indx 

ber (:, k) = (b2 (:, k) -bn (k)) /bn (k) ; 

end 

snr= 0:2:20; 
figure 

subplot(311),plot(snr,100*mean(abs(aer')),'k- 
', snr,100*std(abs(aer')),'k:'),axis([0 20 0 30]) 
title('Normalized Error for a') 
ll=legend('Mean','Std Dev'); 
leg=findobj(11,'type','text'); 

subplot(312),plot(snr,100*meanl(abs(ber')),'k- 
' , snr, 100*stdl(abs(ber')),'k:') ,axis([0 20 0 60]) 
title('Normalized Error for b'), ylabel('Error in %') 

subplot(313),plot(snr,100*mean(abs(xOer')),'k- 
', snr,100*std(abs(xOer')),'k:'),axis([0 20 0 30]) 
title('Normalized Error for xO'),xlabel('SNR') 
toe 
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